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.
asmbgc.F90 in branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_BGC_DA/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10667

Last change on this file since 10667 was 10667, checked in by dford, 6 years ago

Reference FABM indices to jp_fabm_m1.

File size: 144.7 KB
Line 
1MODULE asmbgc
2   !!===========================================================================
3   !!                       ***  MODULE asmbgc  ***
4   !! Routines and declarations for biogeochemical assimilation
5   !!===========================================================================
6   !! History :  3.6  ! 2018-02 (D. Ford)  Adapted from existing modules
7   !!---------------------------------------------------------------------------
8   !! 'key_asminc'          : assimilation increment interface
9   !! 'key_top'             : passive tracers
10   !! 'key_hadocc'          : HadOCC model
11   !! 'key_medusa'          : MEDUSA model
12   !! 'key_roam'            : MEDUSA extras for carbonate cycle
13   !! 'key_fabm'            : ERSEM model coupled via FABM
14   !! 'key_karaml'          : Kara mixed layer depth
15   !!---------------------------------------------------------------------------
16   !! asm_bgc_check_options : check bgc assimilation options
17   !! asm_bgc_init_incs     : initialise bgc increments
18   !! asm_bgc_init_bkg      : initialise bgc background
19   !! asm_bgc_bal_wri       : write out bgc balancing increments
20   !! asm_bgc_bkg_wri       : write out bgc background
21   !! phyto2d_asm_inc       : apply the ocean colour increments
22   !! phyto3d_asm_inc       : apply the 3D phytoplankton increments
23   !! pco2_asm_inc          : apply the pCO2/fCO2 increments
24   !! ph_asm_inc            : apply the pH increments
25   !! bgc3d_asm_inc         : apply the generic 3D BGC increments
26   !!---------------------------------------------------------------------------
27   USE par_kind, ONLY:      & ! kind parameters
28      & wp
29   USE par_oce, ONLY:       & ! domain array sizes
30      & jpi,                &
31      & jpj,                &
32      & jpk
33   USE iom                    ! i/o
34   USE oce, ONLY:           & ! active tracer variables
35      & tsn
36   USE zdfmxl, ONLY :       & ! mixed layer depth
37#if defined key_karaml
38      & hmld_kara,          &
39      & ln_kara,            &
40#endif   
41      & hmld,               &
42      & hmld_tref,          &
43      & hmlp,               &
44      & hmlpt 
45   USE asmpar, ONLY:        & ! assimilation parameters
46      & c_asmbkg,           &
47      & c_asmbal,           &
48      & nitbkg_r,           &
49      & nitavgbkg_r,        &
50      & nitdin_r,           &
51      & nitiaustr_r,        &
52      & nitiaufin_r
53   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
54      & asm_phyto2d_bal_medusa
55   USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
56      & asm_phyto2d_bal_hadocc
57   USE asmphyto2dbal_ersem,  ONLY: & ! phyto2d balancing for ERSEM
58      & asm_phyto2d_bal_ersem
59   USE asmpco2bal, ONLY:    & ! pCO2 balancing
60      & asm_pco2_bal
61#if defined key_top
62   USE trc, ONLY:           & ! passive tracer variables
63      & trn,                &
64      & trb,                &
65      & nittrc000
66   USE par_trc, ONLY:       & ! passive tracer parameters
67      & jptra
68#endif
69#if defined key_medusa
70   USE sms_medusa             ! MEDUSA parameters
71   USE par_medusa             ! MEDUSA parameters
72   USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system
73      & mocsy_interface
74   USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion
75      & f2pCO2
76   USE sms_medusa, ONLY:    & ! MEDUSA diagnostic variables
77      & pgrow_avg,          &
78      & ploss_avg,          &
79      & phyt_avg,           &
80      & mld_max
81#elif defined key_hadocc
82   USE par_hadocc             ! HadOCC parameters
83   USE had_bgc_const          ! HadOCC parameters
84   USE trc, ONLY:           & ! HadOCC diagnostic variables
85      & pgrow_avg,          &
86      & ploss_avg,          &
87      & phyt_avg,           &
88      & mld_max,            &
89      & HADOCC_CHL
90   USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
91      & cchl_p
92#elif defined key_fabm
93   USE par_fabm               ! FABM-ERSEM parameters
94   USE trc, ONLY:           & ! FABM-ERSEM diagnostic variables
95      & pgrow_avg,          &
96      & ploss_avg,          &
97      & phyt_avg,           &
98      & mld_max
99#endif
100
101   IMPLICIT NONE
102   PRIVATE                   
103
104   PUBLIC  asm_bgc_check_options  ! called by asm_inc_init in asminc.F90
105   PUBLIC  asm_bgc_init_incs      ! called by asm_inc_init in asminc.F90
106   PRIVATE asm_bgc_read_incs_2d   ! called by asm_bgc_init_incs
107   PRIVATE asm_bgc_read_incs_3d   ! called by asm_bgc_init_incs
108   PUBLIC  asm_bgc_init_bkg       ! called by asm_inc_init in asminc.F90
109   PUBLIC  asm_bgc_bal_wri        ! called by nemo_gcm in nemogcm.F90
110   PUBLIC  asm_bgc_bkg_alloc      ! called by asm_bkg_wri in asmbkg.F90
111   PUBLIC  asm_bgc_bkg_tavg       ! called by asm_bkg_wri in asmbkg.F90
112   PUBLIC  asm_bgc_bkg_wri        ! called by asm_bkg_wri in asmbkg.F90
113   PUBLIC  phyto2d_asm_inc        ! called by bgc_asm_inc in asminc.F90
114   PUBLIC  phyto3d_asm_inc        ! called by bgc_asm_inc in asminc.F90
115   PUBLIC  pco2_asm_inc           ! called by bgc_asm_inc in asminc.F90
116   PUBLIC  ph_asm_inc             ! called by bgc_asm_inc in asminc.F90
117   PUBLIC  bgc3d_asm_inc          ! called by bgc_asm_inc in asminc.F90
118
119   LOGICAL, PUBLIC :: ln_balwri      = .FALSE. !: No output of balancing incs
120   LOGICAL, PUBLIC :: ln_phytobal    = .FALSE. !: No phytoplankton balancing
121   LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total      log10(chlorophyll) increment
122   LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom     log10(chlorophyll) increment
123   LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
124   LOGICAL, PUBLIC :: ln_slchlnaninc = .FALSE. !: No surface nanophyto  log10(chlorophyll) increment
125   LOGICAL, PUBLIC :: ln_slchlpicinc = .FALSE. !: No surface picophyto  log10(chlorophyll) increment
126   LOGICAL, PUBLIC :: ln_slchldininc = .FALSE. !: No surface dinoflag   log10(chlorophyll) increment
127   LOGICAL, PUBLIC :: ln_schltotinc  = .FALSE. !: No surface total      chlorophyll        increment
128   LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total      log10(phyto C)     increment
129   LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom     log10(phyto C)     increment
130   LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C)     increment
131   LOGICAL, PUBLIC :: ln_spco2inc    = .FALSE. !: No surface pCO2                          increment
132   LOGICAL, PUBLIC :: ln_sfco2inc    = .FALSE. !: No surface fCO2                          increment
133   LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total      log10(chlorophyll) increment
134   LOGICAL, PUBLIC :: ln_pchltotinc  = .FALSE. !: No profile total      chlorophyll        increment
135   LOGICAL, PUBLIC :: ln_pno3inc     = .FALSE. !: No profile nitrate                       increment
136   LOGICAL, PUBLIC :: ln_psi4inc     = .FALSE. !: No profile silicate                      increment
137   LOGICAL, PUBLIC :: ln_ppo4inc     = .FALSE. !: No profile phosphate                     increment
138   LOGICAL, PUBLIC :: ln_pdicinc     = .FALSE. !: No profile dissolved inorganic carbon    increment
139   LOGICAL, PUBLIC :: ln_palkinc     = .FALSE. !: No profile alkalinity                    increment
140   LOGICAL, PUBLIC :: ln_pphinc      = .FALSE. !: No profile pH                            increment
141   LOGICAL, PUBLIC :: ln_po2inc      = .FALSE. !: No profile oxygen                        increment
142
143   INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation
144                                         !  1) hmld      - Turbocline/mixing depth
145                                         !                 [W points]
146                                         !  2) hmlp      - Density criterion
147                                         !                 (0.01 kg/m^3 change from 10m)
148                                         !                 [W points]
149                                         !  3) hmld_kara - Kara MLD
150                                         !                 [Interpolated]
151                                         !  4) hmld_tref - Temperature criterion
152                                         !                 (0.2 K change from surface)
153                                         !                 [T points]
154                                         !  5) hmlpt     - Density criterion
155                                         !                 (0.01 kg/m^3 change from 10m)
156                                         !                 [T points]
157
158   REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
159                                             ! <= 0 no maximum applied (switch turned off)
160                                             !  > 0 absolute chl inc capped at this value
161
162   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchltot_bkginc ! slchltot inc
163   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldia_bkginc ! slchldia inc
164   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
165   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnan_bkginc ! slchlnan inc
166   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlpic_bkginc ! slchlpic inc
167   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldin_bkginc ! slchldin inc
168   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: schltot_bkginc  ! schltot inc
169   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphytot_bkginc ! slphytot inc
170   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphydia_bkginc ! slphydia inc
171   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphynon_bkginc ! slphynon inc
172   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: sfco2_bkginc    ! sfco2 inc
173   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: spco2_bkginc    ! spco2 inc
174   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: plchltot_bkginc ! plchltot inc
175   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pchltot_bkginc  ! pchltot inc
176   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pno3_bkginc     ! pno3 inc
177   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: psi4_bkginc     ! psi4 inc
178   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: ppo4_bkginc     ! ppo4 inc
179   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pdic_bkginc     ! pdic inc
180   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: palk_bkginc     ! palk inc
181   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pph_bkginc      ! pph inc
182   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: po2_bkginc      ! po2 inc
183   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto2d_balinc  ! Balancing incs from ocean colour
184   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto3d_balinc  ! Balancing incs from p(l)chltot
185   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc     ! Balancing incs from spco2/sfco2
186   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc       ! Balancing incs from pph
187
188   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg   ! Background phyto growth
189   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg   ! Background phyto loss
190   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg    ! Background phyto conc
191   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg     ! Background max MLD
192   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg      ! Background tracer state
193   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: chl_bkg         ! Background chl
194   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: cchl_p_bkg      ! Background c:chl
195#if defined key_fabm
196   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: totalk_bkg      ! Background total alkalinity
197#endif
198
199   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: pgrow_avg_tavg  ! Time average of pgrow_avg
200   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: ploss_avg_tavg  ! Time average of ploss_avg
201   REAL(wp), DIMENSION(:,:),     SAVE, ALLOCATABLE :: phyt_avg_tavg   ! Time average of phyt_avg
202   REAL(wp), DIMENSION(:,:,:,:), SAVE, ALLOCATABLE :: trn_tavg        ! Time average of trn
203#if defined key_hadocc
204   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: HADOCC_CHL_tavg ! Time average of HADOCC_CHL
205   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: cchl_p_tavg     ! Time average of cchl_p
206#elif defined key_fabm
207   REAL(wp), DIMENSION(:,:,:),   SAVE, ALLOCATABLE :: totalk_tavg     ! Time average of O3_TA
208#endif
209
210#  include "domzgr_substitute.h90"
211
212CONTAINS
213
214   SUBROUTINE asm_bgc_check_options
215      !!------------------------------------------------------------------------
216      !!                    ***  ROUTINE asm_bgc_check_options  ***
217      !!
218      !! ** Purpose :   check validity of biogeochemical assimilation options
219      !!
220      !! ** Method  :   check settings of logicals
221      !!
222      !! ** Action  :   call ctl_stop/ctl_warn if required
223      !!
224      !! References :   asm_inc_init
225      !!------------------------------------------------------------------------
226
227#if ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa && ! defined key_fabm )
228      CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
229         &           ' but no compatible biogeochemical model is available' )
230#endif
231
232#if defined key_hadocc
233      IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. &
234         & ln_slchldininc .OR. ln_slchlpicinc .OR. ln_slphydiainc .OR. &
235         & ln_slphynoninc .OR. ln_psi4inc     .OR. ln_pphinc      .OR. &
236         & ln_po2inc      .OR. ln_ppo4inc ) THEN
237         CALL ctl_stop( ' Cannot assimilate PFTs, Si4, PO4, pH or O2 into HadOCC' )
238      ENDIF
239#endif
240
241#if defined key_medusa
242      IF ( .NOT. ln_foam_medusa ) THEN
243         CALL ctl_stop( ' MEDUSA data assimilation options not turned on: set ln_foam_medusa' )
244      ENDIF
245
246      IF ( ln_ppo4inc ) THEN
247         CALL ctl_stop( ' Cannot assimilate PO4 into MEDUSA' )
248      ENDIF
249#endif
250
251#if defined key_fabm
252      IF ( ln_pphinc ) THEN
253         CALL ctl_stop( ' Cannot currently assimilate pH into FABM-ERSEM' )
254      ENDIF
255#endif
256
257      IF ( ( ln_phytobal ).AND.                                       &
258         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
259         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
260         & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
261         & ( .NOT. ln_schltotinc  ).AND.( .NOT. ln_slphytotinc ).AND. &
262         & ( .NOT. ln_slphydiainc ).AND.( .NOT. ln_slphynoninc ) ) THEN
263         CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
264            &           ' if not assimilating ocean colour,',                   &
265            &           ' so ln_phytobal will be set to .false.')
266         ln_phytobal = .FALSE.
267      ENDIF
268
269      IF ( ( ln_balwri ).AND.                                         &
270         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
271         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_slchlnaninc ).AND. &
272         & ( .NOT. ln_slchlpicinc ).AND.( .NOT. ln_slchldininc ).AND. &
273         & ( .NOT. ln_schltotinc  ).AND.                              &
274         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
275         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
276         & ( .NOT. ln_pchltotinc  ).AND.( .NOT. ln_pphinc      ).AND. &
277         & ( .NOT. ln_spco2inc    ).AND.( .NOT. ln_sfco2inc    ) ) THEN
278         CALL ctl_warn( ' Balancing increments not being calculated', &
279            &           ' so ln_balwri will be set to .false.')
280         ln_balwri = .FALSE.
281      ENDIF
282
283      IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN
284         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
285      ENDIF
286
287      IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN
288         CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' )
289      ENDIF
290
291      IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN
292         CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' )
293      ENDIF
294
295      IF ( ( ln_slchltotinc .OR. ln_schltotinc  ) .AND. &
296         & ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slchlnaninc .OR. &
297         &   ln_slchlpicinc .OR. ln_slchldininc ) ) THEN
298         CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
299      ENDIF
300
301      IF ( ln_phytobal .AND.                                      &
302         & ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. &
303#if defined key_fabm
304         &   ( (  ( ln_slchldiainc .OR. ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc ) .AND. &
305         &        ( ( .NOT. ln_slchldiainc ) .OR. ( .NOT. ln_slchlnaninc ) .OR. &
306         &          ( .NOT. ln_slchlpicinc ) .OR. ( .NOT. ln_slchldininc ) ) ) .AND. &
307         &     ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) ) THEN
308#elif defined key_medusa
309         &   ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN
310#endif
311         CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
312            &           ' unless assimilating all model PFTs,')
313      ENDIF
314
315      IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN
316         CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' )
317      ENDIF
318
319      IF ( ln_phytobal .AND.                                      &
320         & ( ( ln_slphynoninc .AND.( .NOT. ln_slphydiainc ) ).OR. &
321         &   ( ln_slphydiainc .AND.( .NOT. ln_slphynoninc ) ) ) ) THEN
322         CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
323            &           ' unless assimilating all model PFTs,')
324      ENDIF
325
326   END SUBROUTINE asm_bgc_check_options
327
328   !!===========================================================================
329   !!===========================================================================
330   !!===========================================================================
331
332   SUBROUTINE asm_bgc_init_incs( knum )
333      !!------------------------------------------------------------------------
334      !!                    ***  ROUTINE asm_bgc_init_incs  ***
335      !!
336      !! ** Purpose :   initialise bgc increments
337      !!
338      !! ** Method  :   allocate arrays and read increments from file
339      !!
340      !! ** Action  :   allocate arrays and read increments from file
341      !!
342      !! References :   asm_inc_init
343      !!------------------------------------------------------------------------
344      !!
345      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
346      !!
347      !!------------------------------------------------------------------------
348
349      ! Allocate and read increments
350     
351      IF ( ln_slchltotinc ) THEN
352         ALLOCATE( slchltot_bkginc(jpi,jpj) )
353         CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
354      ENDIF
355     
356      IF ( ln_slchldiainc ) THEN
357         ALLOCATE( slchldia_bkginc(jpi,jpj) )
358         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
359      ENDIF
360     
361      IF ( ln_slchlnoninc ) THEN
362         ALLOCATE( slchlnon_bkginc(jpi,jpj) )
363         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
364      ENDIF
365     
366      IF ( ln_slchlnaninc ) THEN
367         ALLOCATE( slchlnan_bkginc(jpi,jpj) )
368         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnan', slchlnan_bkginc )
369      ENDIF
370     
371      IF ( ln_slchlpicinc ) THEN
372         ALLOCATE( slchlpic_bkginc(jpi,jpj) )
373         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlpic', slchlpic_bkginc )
374      ENDIF
375     
376      IF ( ln_slchldininc ) THEN
377         ALLOCATE( slchldin_bkginc(jpi,jpj) )
378         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldin', slchldin_bkginc )
379      ENDIF
380     
381      IF ( ln_schltotinc ) THEN
382         ALLOCATE( schltot_bkginc(jpi,jpj) )
383         CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
384      ENDIF
385     
386      IF ( ln_slphytotinc ) THEN
387         ALLOCATE( slphytot_bkginc(jpi,jpj) )
388         CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
389      ENDIF
390     
391      IF ( ln_slphydiainc ) THEN
392         ALLOCATE( slphydia_bkginc(jpi,jpj) )
393         CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
394      ENDIF
395     
396      IF ( ln_slphynoninc ) THEN
397         ALLOCATE( slphynon_bkginc(jpi,jpj) )
398         CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc )
399      ENDIF
400
401      IF ( ln_sfco2inc ) THEN
402         ALLOCATE( sfco2_bkginc(jpi,jpj) )
403         CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc )
404      ENDIF
405
406      IF ( ln_spco2inc ) THEN
407         ALLOCATE( spco2_bkginc(jpi,jpj) )
408         CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
409      ENDIF
410     
411      IF ( ln_plchltotinc ) THEN
412         ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
413         CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
414      ENDIF
415     
416      IF ( ln_pchltotinc ) THEN
417         ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
418         CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
419      ENDIF
420     
421      IF ( ln_pno3inc ) THEN
422         ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
423         CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
424      ENDIF
425     
426      IF ( ln_psi4inc ) THEN
427         ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
428         CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
429      ENDIF
430     
431      IF ( ln_ppo4inc ) THEN
432         ALLOCATE( ppo4_bkginc(jpi,jpj,jpk) )
433         CALL asm_bgc_read_incs_3d( knum, 'bckinppo4', ppo4_bkginc )
434      ENDIF
435     
436      IF ( ln_pdicinc ) THEN
437         ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
438         CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
439      ENDIF
440     
441      IF ( ln_palkinc ) THEN
442         ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
443         CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
444      ENDIF
445     
446      IF ( ln_pphinc ) THEN
447         ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
448         CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
449      ENDIF
450     
451      IF ( ln_po2inc ) THEN
452         ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
453         CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc )
454      ENDIF
455
456      ! Allocate balancing increments
457     
458      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
459         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
460         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
461         & ln_slphynoninc ) THEN
462#if defined key_top
463         ALLOCATE( phyto2d_balinc(jpi,jpj,jpk,jptra) )
464         phyto2d_balinc(:,:,:,:) = 0.0
465#else
466         CALL ctl_stop( ' key_top must be set for balancing increments' )
467#endif
468      ENDIF
469     
470      IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
471#if defined key_top
472         ALLOCATE( phyto3d_balinc(jpi,jpj,jpk,jptra) )
473         phyto3d_balinc(:,:,:,:) = 0.0
474#else
475         CALL ctl_stop( ' key_top must be set for balancing increments' )
476#endif
477      ENDIF
478
479      IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
480#if defined key_top
481         ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
482         pco2_balinc(:,:,:,:) = 0.0
483#else
484         CALL ctl_stop( ' key_top must be set for balancing increments' )
485#endif
486      ENDIF
487
488      IF ( ln_pphinc ) THEN
489#if defined key_top
490         ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) )
491         ph_balinc(:,:,:,:) = 0.0
492#else
493         CALL ctl_stop( ' key_top must be set for balancing increments' )
494#endif
495      ENDIF
496
497   END SUBROUTINE asm_bgc_init_incs
498
499   !!===========================================================================
500   !!===========================================================================
501   !!===========================================================================
502
503   SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs )
504      !!------------------------------------------------------------------------
505      !!                    ***  ROUTINE asm_bgc_init_incs  ***
506      !!
507      !! ** Purpose :   read 2d bgc increments
508      !!
509      !! ** Method  :   read increments from file
510      !!
511      !! ** Action  :   read increments from file
512      !!
513      !! References :   asm_inc_init
514      !!------------------------------------------------------------------------
515      !!
516      INTEGER,                      INTENT(in   ) :: knum       ! i/o unit
517      CHARACTER(LEN=*),             INTENT(in   ) :: cd_bgcname ! variable
518      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: p_incs     ! increments
519      !!
520      !!------------------------------------------------------------------------
521
522      ! Initialise
523      p_incs(:,:) = 0.0
524     
525      ! read from file
526      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
527     
528      ! Apply the masks
529      p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
530     
531      ! Set missing increments to 0.0 rather than 1e+20
532      ! to allow for differences in masks
533      WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0
534
535   END SUBROUTINE asm_bgc_read_incs_2d
536
537   !!===========================================================================
538   !!===========================================================================
539   !!===========================================================================
540
541   SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs )
542      !!------------------------------------------------------------------------
543      !!                    ***  ROUTINE asm_bgc_init_incs  ***
544      !!
545      !! ** Purpose :   read 3d bgc increments
546      !!
547      !! ** Method  :   read increments from file
548      !!
549      !! ** Action  :   read increments from file
550      !!
551      !! References :   asm_inc_init
552      !!------------------------------------------------------------------------
553      !!
554      INTEGER,                          INTENT(in   ) :: knum       ! i/o unit
555      CHARACTER(LEN=*),                 INTENT(in   ) :: cd_bgcname ! variable
556      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: p_incs     ! increments
557      !!
558      !!------------------------------------------------------------------------
559
560      ! Initialise
561      p_incs(:,:,:) = 0.0
562     
563      ! read from file
564      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
565     
566      ! Apply the masks
567      p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
568     
569      ! Set missing increments to 0.0 rather than 1e+20
570      ! to allow for differences in masks
571      WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0
572
573   END SUBROUTINE asm_bgc_read_incs_3d
574
575   !!===========================================================================
576   !!===========================================================================
577   !!===========================================================================
578
579   SUBROUTINE asm_bgc_init_bkg
580      !!------------------------------------------------------------------------
581      !!                    ***  ROUTINE asm_bgc_init_bkg  ***
582      !!
583      !! ** Purpose :   initialise bgc background
584      !!
585      !! ** Method  :   allocate arrays and read background from file
586      !!
587      !! ** Action  :   allocate arrays and read background from file
588      !!
589      !! References :   asm_inc_init
590      !!------------------------------------------------------------------------
591      !!
592      INTEGER :: inum   ! i/o unit of background file
593      INTEGER :: jt     ! loop counter
594      !!
595      !!------------------------------------------------------------------------
596
597#if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm )
598      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
599         & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
600         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
601         & ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN
602
603         ALLOCATE( pgrow_avg_bkg(jpi,jpj)        )
604         ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
605         ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
606         ALLOCATE( mld_max_bkg(jpi,jpj)          )
607         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
608         pgrow_avg_bkg(:,:)  = 0.0
609         ploss_avg_bkg(:,:)  = 0.0
610         phyt_avg_bkg(:,:)   = 0.0
611         mld_max_bkg(:,:)    = 0.0
612         tracer_bkg(:,:,:,:) = 0.0
613
614#if defined key_hadocc
615         ALLOCATE( chl_bkg(jpi,jpj,jpk)    )
616         ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) )
617         chl_bkg(:,:,:)    = 0.0
618         cchl_p_bkg(:,:,:) = 0.0
619#elif defined key_fabm
620         ALLOCATE( totalk_bkg(jpi,jpj,jpk) )
621         totalk_bkg(:,:,:) = 0.0
622#endif
623         
624         !--------------------------------------------------------------------
625         ! Read background variables for phytoplankton assimilation
626         ! Some only required if performing balancing
627         !--------------------------------------------------------------------
628
629         CALL iom_open( c_asmbkg, inum )
630
631#if defined key_hadocc
632         CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl',  chl_bkg    )
633         CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
634         chl_bkg(:,:,:)    = chl_bkg(:,:,:)    * tmask(:,:,:)
635         cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:)
636#elif defined key_medusa
637         CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
638         CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
639         CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
640         CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
641         CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
642#elif defined key_fabm
643         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl1', tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl1) )
644         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl2', tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl2) )
645         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl3', tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl3) )
646         CALL iom_get( inum, jpdom_autoglo, 'ersem_chl4', tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl4) )
647         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p1c)  )
648         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p1n)  )
649         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1p',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p1p)  )
650         CALL iom_get( inum, jpdom_autoglo, 'ersem_p1s',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p1s)  )
651         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p2c)  )
652         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p2n)  )
653         CALL iom_get( inum, jpdom_autoglo, 'ersem_p2p',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p2p)  )
654         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p3c)  )
655         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p3n)  )
656         CALL iom_get( inum, jpdom_autoglo, 'ersem_p3p',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p3p)  )
657         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p4c)  )
658         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p4n)  )
659         CALL iom_get( inum, jpdom_autoglo, 'ersem_p4p',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_p4p)  )
660#endif
661         
662         IF ( ln_phytobal ) THEN
663
664            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
665            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
666            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  )
667            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
668            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
669            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
670            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1)
671            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1)
672
673#if defined key_hadocc
674            CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
675            CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
676            CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
677            CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
678            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
679            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
680#elif defined key_medusa
681            CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
682            CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
683            CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
684            CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
685            CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
686            CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
687            CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
688            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
689            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
690            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
691#elif defined key_fabm
692            CALL iom_get( inum, jpdom_autoglo, 'ersem_z4c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z4c)  )
693            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z5c)  )
694            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z5n)  )
695            CALL iom_get( inum, jpdom_autoglo, 'ersem_z5p',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z5p)  )
696            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z6c)  )
697            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z6n)  )
698            CALL iom_get( inum, jpdom_autoglo, 'ersem_z6p',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_z6p)  )
699            CALL iom_get( inum, jpdom_autoglo, 'ersem_n1p',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n1p)  )
700            CALL iom_get( inum, jpdom_autoglo, 'ersem_n3n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n3n)  )
701            CALL iom_get( inum, jpdom_autoglo, 'ersem_n4n',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n4n)  )
702            CALL iom_get( inum, jpdom_autoglo, 'ersem_n5s',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_n5s)  )
703            CALL iom_get( inum, jpdom_autoglo, 'ersem_o2o',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o2o)  )
704            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
705            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
706            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              )
707            totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
708#endif
709         ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
710#if defined key_hadocc
711            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
712            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
713#elif defined key_medusa
714            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
715            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
716#elif defined key_fabm
717            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
718            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
719            CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              )
720            totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
721#endif
722            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
723            mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
724         ENDIF
725
726         CALL iom_close( inum )
727         
728         DO jt = 1, jptra
729            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
730         END DO
731     
732      ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
733
734         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
735         ALLOCATE( mld_max_bkg(jpi,jpj)          )
736         tracer_bkg(:,:,:,:) = 0.0
737         mld_max_bkg(:,:)    = 0.0
738
739         CALL iom_open( c_asmbkg, inum )
740         
741#if defined key_hadocc
742         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
743         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
744#elif defined key_medusa
745         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
746         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
747#elif defined key_fabm
748         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3c',  tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
749         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ba', tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
750         CALL iom_get( inum, jpdom_autoglo, 'ersem_o3ta', totalk_bkg(:,:,:)              )
751         totalk_bkg(:,:,:) = totalk_bkg(:,:,:) * tmask(:,:,:)
752#endif
753         CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
754
755         CALL iom_close( inum )
756         
757         DO jt = 1, jptra
758            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
759         END DO
760         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
761     
762      ENDIF
763#else
764      CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' )
765#endif
766
767   END SUBROUTINE asm_bgc_init_bkg
768
769   !!===========================================================================
770   !!===========================================================================
771   !!===========================================================================
772
773   SUBROUTINE asm_bgc_bal_wri( kt )
774      !!------------------------------------------------------------------------
775      !!
776      !!                  ***  ROUTINE asm_bgc_bal_wri ***
777      !!
778      !! ** Purpose : Write to file the balancing increments
779      !!              calculated for biogeochemistry
780      !!
781      !! ** Method  : Write to file the balancing increments
782      !!              calculated for biogeochemistry
783      !!
784      !! ** Action  :
785      !!                   
786      !! References :
787      !!
788      !! History    :
789      !!        ! 2014-08 (D. Ford)  Initial version, based on asm_bkg_wri
790      !!------------------------------------------------------------------------
791      !! * Arguments
792      INTEGER, INTENT( IN ) :: kt        ! Current time-step
793      !! * Local declarations
794      CHARACTER(LEN=50) :: cl_asmbal     ! Filename (with extension)
795      LOGICAL           :: llok          ! Check if file exists
796      INTEGER           :: inum          ! File unit number
797      REAL(wp)          :: zdate         ! Date
798      !!------------------------------------------------------------------------
799     
800      ! Set things up
801      zdate = REAL( ndastp )
802      WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
803
804      ! Check if file exists
805      INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
806      IF( .NOT. llok ) THEN
807         IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
808            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
809
810         ! Define the output file       
811         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
812
813         ! Write the information
814         CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   )
815
816         IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
817            & ln_slchlnaninc .OR. ln_slchlpicinc .OR. ln_slchldininc .OR. &
818            & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
819            & ln_slphynoninc ) THEN
820#if defined key_medusa
821            CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) )
822            CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) )
823            CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) )
824            CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) )
825            CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) )
826            IF ( ln_phytobal ) THEN
827               CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) )
828               CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) )
829               CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) )
830               CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) )
831               CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) )
832               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) )
833               CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) )
834               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) )
835               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) )
836               CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) )
837            ENDIF
838#elif defined key_hadocc
839            CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
840            IF ( ln_phytobal ) THEN
841               CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
842               CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
843               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
844               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
845               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
846            ENDIF
847#elif defined key_fabm
848            CALL iom_rstput( kt, kt, inum, 'phy2d_chl1', phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl1) )
849            CALL iom_rstput( kt, kt, inum, 'phy2d_chl2', phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl2) )
850            CALL iom_rstput( kt, kt, inum, 'phy2d_chl3', phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl3) )
851            CALL iom_rstput( kt, kt, inum, 'phy2d_chl4', phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl4) )
852            CALL iom_rstput( kt, kt, inum, 'phy2d_p1c',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1c)  )
853            CALL iom_rstput( kt, kt, inum, 'phy2d_p1n',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1n)  )
854            CALL iom_rstput( kt, kt, inum, 'phy2d_p1p',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1p)  )
855            CALL iom_rstput( kt, kt, inum, 'phy2d_p1s',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1s)  )
856            CALL iom_rstput( kt, kt, inum, 'phy2d_p2c',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2c)  )
857            CALL iom_rstput( kt, kt, inum, 'phy2d_p2n',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2n)  )
858            CALL iom_rstput( kt, kt, inum, 'phy2d_p2p',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2p)  )
859            CALL iom_rstput( kt, kt, inum, 'phy2d_p3c',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3c)  )
860            CALL iom_rstput( kt, kt, inum, 'phy2d_p3n',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3n)  )
861            CALL iom_rstput( kt, kt, inum, 'phy2d_p3p',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3p)  )
862            CALL iom_rstput( kt, kt, inum, 'phy2d_p4c',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4c)  )
863            CALL iom_rstput( kt, kt, inum, 'phy2d_p4n',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4n)  )
864            CALL iom_rstput( kt, kt, inum, 'phy2d_p4p',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4p)  )
865            IF ( ln_phytobal ) THEN
866               CALL iom_rstput( kt, kt, inum, 'phy2d_z4c',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_z4c)  )
867               CALL iom_rstput( kt, kt, inum, 'phy2d_z5c',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_z5c)  )
868               CALL iom_rstput( kt, kt, inum, 'phy2d_z5n',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_z5n)  )
869               CALL iom_rstput( kt, kt, inum, 'phy2d_z5p',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_z5p)  )
870               CALL iom_rstput( kt, kt, inum, 'phy2d_z6c',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_z6c)  )
871               CALL iom_rstput( kt, kt, inum, 'phy2d_z6n',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_z6n)  )
872               CALL iom_rstput( kt, kt, inum, 'phy2d_z6p',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_z6p)  )
873               CALL iom_rstput( kt, kt, inum, 'phy2d_n1p',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_n1p)  )
874               CALL iom_rstput( kt, kt, inum, 'phy2d_n3n',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_n3n)  )
875               CALL iom_rstput( kt, kt, inum, 'phy2d_n4n',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_n4n)  )
876               CALL iom_rstput( kt, kt, inum, 'phy2d_n5s',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_n5s)  )
877               CALL iom_rstput( kt, kt, inum, 'phy2d_o2o',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_o2o)  )
878               CALL iom_rstput( kt, kt, inum, 'phy2d_o3c',   phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
879               CALL iom_rstput( kt, kt, inum, 'phy2d_o3ba',  phyto2d_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
880            ENDIF
881#endif
882         ENDIF
883
884         IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
885#if defined key_medusa
886            CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) )
887            CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) )
888            CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) )
889            CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
890            CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
891#elif defined key_hadocc
892            CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
893#elif defined key_fabm
894            CALL iom_rstput( kt, kt, inum, 'phy3d_chl1', phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl1) )
895            CALL iom_rstput( kt, kt, inum, 'phy3d_chl2', phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl2) )
896            CALL iom_rstput( kt, kt, inum, 'phy3d_chl3', phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl3) )
897            CALL iom_rstput( kt, kt, inum, 'phy3d_chl4', phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_chl4) )
898            CALL iom_rstput( kt, kt, inum, 'phy3d_p1c',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1c)  )
899            CALL iom_rstput( kt, kt, inum, 'phy3d_p1n',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1n)  )
900            CALL iom_rstput( kt, kt, inum, 'phy3d_p1p',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1p)  )
901            CALL iom_rstput( kt, kt, inum, 'phy3d_p1s',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p1s)  )
902            CALL iom_rstput( kt, kt, inum, 'phy3d_p2c',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2c)  )
903            CALL iom_rstput( kt, kt, inum, 'phy3d_p2n',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2n)  )
904            CALL iom_rstput( kt, kt, inum, 'phy3d_p2p',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p2p)  )
905            CALL iom_rstput( kt, kt, inum, 'phy3d_p3c',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3c)  )
906            CALL iom_rstput( kt, kt, inum, 'phy3d_p3n',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3n)  )
907            CALL iom_rstput( kt, kt, inum, 'phy3d_p3p',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p3p)  )
908            CALL iom_rstput( kt, kt, inum, 'phy3d_p4c',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4c)  )
909            CALL iom_rstput( kt, kt, inum, 'phy3d_p4n',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4n)  )
910            CALL iom_rstput( kt, kt, inum, 'phy3d_p4p',  phyto3d_balinc(:,:,:,jp_fabm_m1+jp_fabm_p4p)  )
911#endif
912         ENDIF
913
914         IF ( ln_spco2inc ) THEN
915#if defined key_medusa
916            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) )
917            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) )
918#elif defined key_hadocc
919            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
920            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
921#elif defined key_fabm
922            CALL iom_rstput( kt, kt, inum, 'pco2_o3c',  pco2_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
923            CALL iom_rstput( kt, kt, inum, 'pco2_o3ba', pco2_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
924#endif
925         ELSE IF ( ln_sfco2inc ) THEN
926#if defined key_medusa
927            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) )
928            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) )
929#elif defined key_hadocc
930            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
931            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
932#elif defined key_fabm
933            CALL iom_rstput( kt, kt, inum, 'fco2_o3c',  pco2_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
934            CALL iom_rstput( kt, kt, inum, 'fco2_o3ba', pco2_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
935#endif
936         ENDIF
937
938         IF ( ln_pphinc ) THEN
939#if defined key_medusa
940            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) )
941            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) )
942#elif defined key_hadocc
943            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
944            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
945#elif defined key_fabm
946            CALL iom_rstput( kt, kt, inum, 'ph_o3c',  ph_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
947            CALL iom_rstput( kt, kt, inum, 'ph_o3ba', ph_balinc(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
948#endif
949         ENDIF
950
951         CALL iom_close( inum )
952      ELSE
953         CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
954            &           ' Therefore not writing out balancing increments at this timestep', &
955            &           ' Something has probably gone wrong somewhere' )
956         IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
957            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
958      ENDIF
959                                 
960   END SUBROUTINE asm_bgc_bal_wri
961
962   !!===========================================================================
963   !!===========================================================================
964   !!===========================================================================
965
966   SUBROUTINE asm_bgc_bkg_alloc
967      !!------------------------------------------------------------------------
968      !!                    ***  ROUTINE asm_bgc_bkg_alloc  ***
969      !!
970      !! ** Purpose :   allocate time-average arrays for background
971      !!
972      !! ** Method  :   allocate time-average arrays for background
973      !!
974      !! ** Action  :   allocate time-average arrays for background
975      !!
976      !! References :   asm_bkg_wri
977      !!------------------------------------------------------------------------
978      !!
979      INTEGER :: ierror
980      !!
981      !!------------------------------------------------------------------------
982     
983      ALLOCATE( pgrow_avg_tavg(jpi,jpj), STAT=ierror )
984      IF( ierror > 0 ) THEN
985         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate pgrow_avg_tavg' )
986      ENDIF
987      pgrow_avg_tavg(:,:) = 0.0
988     
989      ALLOCATE( ploss_avg_tavg(jpi,jpj), STAT=ierror )
990      IF( ierror > 0 ) THEN
991         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate ploss_avg_tavg' )
992      ENDIF
993      ploss_avg_tavg(:,:) = 0.0
994     
995      ALLOCATE( phyt_avg_tavg(jpi,jpj), STAT=ierror )
996      IF( ierror > 0 ) THEN
997         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate phyt_avg_tavg' )
998      ENDIF
999      phyt_avg_tavg(:,:) = 0.0
1000     
1001      ALLOCATE( trn_tavg(jpi,jpj,jpk,jptra), STAT=ierror )
1002      IF( ierror > 0 ) THEN
1003         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate trn_tavg' )
1004      ENDIF
1005      trn_tavg(:,:,:,:) = 0.0
1006
1007#if defined key_hadocc
1008      ALLOCATE( HADOCC_CHL_tavg(jpi,jpj,jpk), STAT=ierror )
1009      IF( ierror > 0 ) THEN
1010         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate HADOCC_CHL_tavg' )
1011      ENDIF
1012      HADOCC_CHL_tavg(:,:,:) = 0.0
1013
1014      ALLOCATE( cchl_p_tavg(jpi,jpj,jpk), STAT=ierror )
1015      IF( ierror > 0 ) THEN
1016         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate cchl_p_tavg' )
1017      ENDIF
1018      cchl_p_tavg(:,:,:) = 0.0
1019
1020#elif defined key_fabm
1021      ALLOCATE( totalk_tavg(jpi,jpj,jpk), STAT=ierror )
1022      IF( ierror > 0 ) THEN
1023         CALL ctl_stop( 'asm_bgc_bkg_alloc: unable to allocate totalk_tavg' )
1024      ENDIF
1025      totalk_tavg(:,:,:) = 0.0
1026#endif
1027                                 
1028   END SUBROUTINE asm_bgc_bkg_alloc
1029
1030   !!===========================================================================
1031   !!===========================================================================
1032   !!===========================================================================
1033
1034   SUBROUTINE asm_bgc_bkg_tavg(kt, pnumtimes_tavg)
1035      !!------------------------------------------------------------------------
1036      !!                    ***  ROUTINE asm_bgc_bkg_wri  ***
1037      !!
1038      !! ** Purpose :   write out bgc background
1039      !!
1040      !! ** Method  :   write out bgc background
1041      !!
1042      !! ** Action  :   write out bgc background
1043      !!
1044      !! References :   asm_bkg_wri
1045      !!------------------------------------------------------------------------
1046      !!
1047      INTEGER,  INTENT(in   ) :: kt               ! Current time-step
1048      !!
1049      REAL(wp), INTENT(in   ) :: pnumtimes_tavg   ! No of times to average over
1050      !!
1051      !!------------------------------------------------------------------------
1052     
1053      IF (kt == nittrc000) THEN
1054         pgrow_avg_tavg(:,:) = 0.0
1055         ploss_avg_tavg(:,:) = 0.0
1056         phyt_avg_tavg(:,:)  = 0.0
1057         trn_tavg(:,:,:,:)   = 0.0
1058#if defined key_hadocc
1059         HADOCC_CHL_tavg(:,:,:) = 0.0
1060         cchl_p_tavg(:,:,:)     = 0.0
1061#elif defined key_fabm
1062         totalk_tavg(:,:,:)  = 0.0
1063#endif
1064      ENDIF
1065     
1066      pgrow_avg_tavg(:,:) = pgrow_avg_tavg(:,:) + pgrow_avg(:,:) / pnumtimes_tavg
1067      ploss_avg_tavg(:,:) = ploss_avg_tavg(:,:) + ploss_avg(:,:) / pnumtimes_tavg
1068      phyt_avg_tavg(:,:)  = phyt_avg_tavg(:,:)  + phyt_avg(:,:)  / pnumtimes_tavg
1069      trn_tavg(:,:,:,:)   = trn_tavg(:,:,:,:)   + trn(:,:,:,:)   / pnumtimes_tavg
1070#if defined key_hadocc
1071      HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL_tavg(:,:,:) + HADOCC_CHL(:,:,:) / pnumtimes_tavg
1072      cchl_p_tavg(:,:,:)     = cchl_p_tavg(:,:,:)     + cchl_p(:,:,:)     / pnumtimes_tavg
1073#elif defined key_fabm
1074      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) + &
1075         &                  fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta) / pnumtimes_tavg
1076      totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:)
1077#endif
1078                                 
1079   END SUBROUTINE asm_bgc_bkg_tavg
1080
1081   !!===========================================================================
1082   !!===========================================================================
1083   !!===========================================================================
1084
1085   SUBROUTINE asm_bgc_bkg_wri( kt, knum, ld_avgbkg )
1086      !!------------------------------------------------------------------------
1087      !!                    ***  ROUTINE asm_bgc_bkg_wri  ***
1088      !!
1089      !! ** Purpose :   write out bgc background
1090      !!
1091      !! ** Method  :   write out bgc background
1092      !!
1093      !! ** Action  :   write out bgc background
1094      !!
1095      !! References :   asm_bkg_wri
1096      !!------------------------------------------------------------------------
1097      !!
1098      INTEGER, INTENT(in   ) :: kt          ! Current time-step
1099      INTEGER, INTENT(in   ) :: knum        ! i/o unit of increments file
1100      LOGICAL, INTENT(in   ) :: ld_avgbkg   ! Averaged background?
1101      !!
1102      INTEGER                :: nitbgcbkg_r ! Period referenced to nit000
1103      !!
1104      !!------------------------------------------------------------------------
1105     
1106      IF (ld_avgbkg) THEN
1107         nitbgcbkg_r = nitavgbkg_r
1108      ELSE
1109         nitbgcbkg_r = nitbkg_r
1110         pgrow_avg_tavg(:,:) = pgrow_avg(:,:)
1111         ploss_avg_tavg(:,:) = ploss_avg(:,:)
1112         phyt_avg_tavg(:,:)  = phyt_avg(:,:)
1113         trn_tavg(:,:,:,:)   = trn(:,:,:,:)
1114#if defined key_hadocc
1115         HADOCC_CHL_tavg(:,:,:) = HADOCC_CHL(:,:,:)
1116         cchl_p_tavg(:,:,:)     = cchl_p(:,:,:)
1117#elif defined key_fabm
1118         totalk_tavg(:,:,:)  = fabm_get_interior_diagnostic_data(model, jp_fabm_o3ta)
1119         totalk_tavg(:,:,:)  = totalk_tavg(:,:,:) * tmask(:,:,:)
1120#endif
1121      ENDIF
1122
1123#if defined key_hadocc
1124      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg             )
1125      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg             )
1126      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg              )
1127      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max                    )
1128      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_nut'  , trn_tavg(:,:,:,jp_had_nut) )
1129      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_phy'  , trn_tavg(:,:,:,jp_had_phy) )
1130      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_zoo'  , trn_tavg(:,:,:,jp_had_zoo) )
1131      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_pdn'  , trn_tavg(:,:,:,jp_had_pdn) )
1132      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_dic'  , trn_tavg(:,:,:,jp_had_dic) )
1133      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_alk'  , trn_tavg(:,:,:,jp_had_alk) )
1134      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,:)     )
1135      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         )
1136#elif defined key_medusa
1137      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg        )
1138      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg        )
1139      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg         )
1140      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max               )
1141      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chn'  , trn_tavg(:,:,:,jpchn) )
1142      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_chd'  , trn_tavg(:,:,:,jpchd) )
1143      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phn'  , trn_tavg(:,:,:,jpphn) )
1144      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_phd'  , trn_tavg(:,:,:,jpphd) )
1145      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_pds'  , trn_tavg(:,:,:,jppds) )
1146      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zmi'  , trn_tavg(:,:,:,jpzmi) )
1147      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_zme'  , trn_tavg(:,:,:,jpzme) )
1148      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_din'  , trn_tavg(:,:,:,jpdin) )
1149      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_sil'  , trn_tavg(:,:,:,jpsil) )
1150      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_fer'  , trn_tavg(:,:,:,jpfer) )
1151      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_det'  , trn_tavg(:,:,:,jpdet) )
1152      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dtc'  , trn_tavg(:,:,:,jpdtc) )
1153      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_dic'  , trn_tavg(:,:,:,jpdic) )
1154      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_alk'  , trn_tavg(:,:,:,jpalk) )
1155      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'medusa_oxy'  , trn_tavg(:,:,:,jpoxy) )
1156#elif defined key_fabm
1157      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'pgrow_avg'   , pgrow_avg_tavg               )
1158      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ploss_avg'   , ploss_avg_tavg               )
1159      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'phyt_avg'    , phyt_avg_tavg                )
1160      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'mld_max'     , mld_max                      )
1161      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl1'  , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_chl1) )
1162      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl2'  , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_chl2) )
1163      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl3'  , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_chl3) )
1164      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_chl4'  , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_chl4) )
1165      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p1c)  )
1166      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p1n)  )
1167      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1p'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p1p)  )
1168      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p1s'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p1s)  )
1169      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p2c)  )
1170      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p2n)  )
1171      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p2p'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p2p)  )
1172      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p3c)  )
1173      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p3n)  )
1174      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p3p'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p3p)  )
1175      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p4c)  )
1176      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p4n)  )
1177      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_p4p'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_p4p)  )
1178      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z4c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_z4c)  )
1179      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_z5c)  )
1180      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_z5n)  )
1181      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z5p'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_z5p)  )
1182      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_z6c)  )
1183      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_z6n)  )
1184      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_z6p'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_z6p)  )
1185      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n1p'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_n1p)  )
1186      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n3n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_n3n)  )
1187      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n4n'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_n4n)  )
1188      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_n5s'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_n5s)  )
1189      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o2o'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_o2o)  )
1190      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3c'   , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_o3c)  )
1191      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ba'  , trn_tavg(:,:,:,jp_fabm_m1+jp_fabm_o3ba) )
1192      CALL iom_rstput( kt, nitbgcbkg_r, knum, 'ersem_o3ta'  , totalk_tavg                  )
1193#endif
1194
1195   END SUBROUTINE asm_bgc_bkg_wri
1196
1197   !!===========================================================================
1198   !!===========================================================================
1199   !!===========================================================================
1200
1201   SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog )
1202      !!------------------------------------------------------------------------
1203      !!                    ***  ROUTINE asm_bgc_init_incs  ***
1204      !!
1205      !! ** Purpose :   convert log increments to non-log
1206      !!
1207      !! ** Method  :   need to account for model background,
1208      !!                cannot simply do 10^log_inc. Need to:
1209      !!                1) Add log_inc to log10(background) to get log10(analysis)
1210      !!                2) Take 10^log10(analysis) to get analysis
1211      !!                3) Subtract background from analysis to get non-log incs
1212      !!
1213      !! ** Action  :   return non-log increments
1214      !!
1215      !! References :   
1216      !!------------------------------------------------------------------------
1217      !!
1218      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pbkg        ! Background
1219      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pinc_log    ! Log incs
1220      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs
1221      !
1222      INTEGER                                     :: ji, jj      ! Loop counters
1223      !!
1224      !!------------------------------------------------------------------------
1225
1226      DO jj = 1, jpj
1227         DO ji = 1, jpi
1228            IF ( pbkg(ji,jj) > 0.0 ) THEN
1229               pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + &
1230                  &                       pinc_log(ji,jj) )      &
1231                  &                 - pbkg(ji,jj)
1232            ELSE
1233               pinc_nonlog(ji,jj) = 0.0
1234            ENDIF
1235         END DO
1236      END DO
1237
1238   END SUBROUTINE asm_bgc_unlog_2d
1239
1240   !!===========================================================================
1241   !!===========================================================================
1242   !!===========================================================================
1243
1244   SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1245      !!------------------------------------------------------------------------
1246      !!                    ***  ROUTINE phyto2d_asm_inc  ***
1247      !!         
1248      !! ** Purpose : Apply the chlorophyll assimilation increments.
1249      !!
1250      !! ** Method  : Calculate increments to state variables using nitrogen
1251      !!              balancing.
1252      !!              Direct initialization or Incremental Analysis Updating.
1253      !!
1254      !! ** Action  :
1255      !!------------------------------------------------------------------------
1256      INTEGER,  INTENT(IN) :: kt        ! Current time step
1257      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1258      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1259      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1260      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1261      !
1262      INTEGER                      :: it          ! Index
1263      REAL(wp)                     :: zincwgt     ! IAU weight for current time step
1264      REAL(wp)                     :: zincper     ! IAU interval in seconds
1265      REAL(wp), DIMENSION(jpi,jpj) :: zmld        ! Mixed layer depth
1266      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs
1267      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
1268      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs
1269      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
1270#if defined key_medusa
1271      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
1272      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
1273      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs
1274      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
1275      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs
1276      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
1277      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
1278      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
1279#elif defined key_fabm
1280      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
1281      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
1282      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldin ! Local chldin incs
1283      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldin ! Local chldin bkg
1284      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnan ! Local chlnan incs
1285      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnan ! Local chlnan bkg
1286      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlpic ! Local chlpic incs
1287      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlpic ! Local chlpic bkg
1288#endif
1289      !!------------------------------------------------------------------------
1290     
1291      IF ( kt <= nit000 ) THEN
1292
1293         ! Un-log any log increments for passing to balancing routines
1294         ! Remember that two sets of non-log increments should not be
1295         ! expected to be in the same ratio as their log equivalents
1296         
1297         ! Total chlorophyll
1298         IF ( ln_slchltotinc ) THEN
1299#if defined key_medusa
1300            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
1301#elif defined key_hadocc
1302            zbkg_chltot(:,:) = chl_bkg(:,:,1)
1303#elif defined key_fabm
1304            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl1) + &
1305               &               tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl2) + &
1306               &               tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl3) + &
1307               &               tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl4)
1308#endif
1309            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
1310         ELSE IF ( ln_schltotinc ) THEN
1311            zinc_chltot(:,:) = schltot_bkginc(:,:)
1312         ELSE
1313            zinc_chltot(:,:) = 0.0
1314         ENDIF
1315
1316#if defined key_medusa || defined key_fabm
1317         ! Diatom chlorophyll
1318         IF ( ln_slchldiainc ) THEN
1319#if defined key_medusa
1320            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
1321#elif defined key_fabm
1322            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl1)
1323#endif
1324            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
1325         ELSE
1326            zinc_chldia(:,:) = 0.0
1327         ENDIF
1328#endif
1329
1330#if defined key_medusa
1331         ! Non-diatom chlorophyll
1332         IF ( ln_slchlnoninc ) THEN
1333            zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
1334            CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon )
1335         ELSE
1336            zinc_chlnon(:,:) = 0.0
1337         ENDIF
1338#endif
1339
1340#if defined key_fabm
1341         ! Nanophytoplankton chlorophyll
1342         IF ( ln_slchlnaninc ) THEN
1343            zbkg_chlnan(:,:) = tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl2)
1344            CALL asm_bgc_unlog_2d( zbkg_chlnan, slchlnan_bkginc, zinc_chlnan )
1345         ELSE
1346            zinc_chlnan(:,:) = 0.0
1347         ENDIF
1348
1349         ! Picophytoplankton chlorophyll
1350         IF ( ln_slchlpicinc ) THEN
1351            zbkg_chlpic(:,:) = tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl3)
1352            CALL asm_bgc_unlog_2d( zbkg_chlpic, slchlpic_bkginc, zinc_chlpic )
1353         ELSE
1354            zinc_chlpic(:,:) = 0.0
1355         ENDIF
1356
1357         ! Dinoflagellate chlorophyll
1358         IF ( ln_slchldininc ) THEN
1359            zbkg_chldin(:,:) = tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_chl4)
1360            CALL asm_bgc_unlog_2d( zbkg_chldin, slchldin_bkginc, zinc_chldin )
1361         ELSE
1362            zinc_chldin(:,:) = 0.0
1363         ENDIF
1364#endif
1365
1366         ! Total phytoplankton carbon
1367         IF ( ln_slphytotinc ) THEN
1368#if defined key_medusa
1369            zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
1370#elif defined key_hadocc
1371            zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
1372#endif
1373            CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot )
1374         ELSE
1375            zinc_phytot(:,:) = 0.0
1376         ENDIF
1377
1378#if defined key_medusa
1379         ! Diatom phytoplankton carbon
1380         IF ( ln_slphydiainc ) THEN
1381            zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
1382            CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia )
1383         ELSE
1384            zinc_phydia(:,:) = 0.0
1385         ENDIF
1386#endif
1387
1388#if defined key_medusa
1389         ! Non-diatom phytoplankton carbon
1390         IF ( ln_slphynoninc ) THEN
1391            zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
1392            CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon )
1393         ELSE
1394            zinc_phynon(:,:) = 0.0
1395         ENDIF
1396#endif
1397
1398         ! Select mixed layer
1399         IF ( ll_asmdin ) THEN
1400#if defined key_top && ( defined key_hadocc || defined key_medusa || defined key_fabm )
1401            CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
1402               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1403            zmld(:,:) = mld_max_bkg(:,:)
1404#else
1405            CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
1406#endif
1407         ELSE
1408            SELECT CASE( mld_choice_bgc )
1409            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
1410               zmld(:,:) = hmld(:,:)
1411            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1412               zmld(:,:) = hmlp(:,:)
1413            CASE ( 3 )                   ! Kara MLD [Interpolated]
1414#if defined key_karaml
1415               IF ( ln_kara ) THEN
1416                  zmld(:,:) = hmld_kara(:,:)
1417               ELSE
1418                  CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1419                     &           ' but ln_kara=.false.' )
1420               ENDIF
1421#else
1422               CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1423                  &           ' but is not defined' )
1424#endif
1425            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
1426               zmld(:,:) = hmld_tref(:,:)
1427            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1428               zmld(:,:) = hmlpt(:,:)
1429            END SELECT
1430         ENDIF
1431
1432         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
1433
1434#if defined key_medusa
1435         CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
1436            &                         zinc_chltot,                         &
1437            &                         ln_slchldiainc,                      &
1438            &                         zinc_chldia,                         &
1439            &                         ln_slchlnoninc,                      &
1440            &                         zinc_chlnon,                         &
1441            &                         ln_slphytotinc,                      &
1442            &                         zinc_phytot,                         &
1443            &                         ln_slphydiainc,                      &
1444            &                         zinc_phydia,                         &
1445            &                         ln_slphynoninc,                      &
1446            &                         zinc_phynon,                         &
1447            &                         zincper,                             &
1448            &                         rn_maxchlinc, ln_phytobal, zmld,     &
1449            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
1450            &                         phyt_avg_bkg, mld_max_bkg,           &
1451            &                         tracer_bkg, phyto2d_balinc )
1452#elif defined key_hadocc
1453         CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
1454            &                         zinc_chltot,                         &
1455            &                         ln_slphytotinc,                      &
1456            &                         zinc_phytot,                         &
1457            &                         zincper,                             &
1458            &                         rn_maxchlinc, ln_phytobal, zmld,     &
1459            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
1460            &                         phyt_avg_bkg, mld_max_bkg,           &
1461            &                         cchl_p_bkg(:,:,1),                   &
1462            &                         tracer_bkg, phyto2d_balinc )
1463#elif defined key_fabm
1464         CALL asm_phyto2d_bal_ersem(  (ln_slchltotinc .OR. ln_schltotinc), &
1465            &                         zinc_chltot,                         &
1466            &                         ln_slchldiainc,                      &
1467            &                         zinc_chldia,                         &
1468            &                         ln_slchlnaninc,                      &
1469            &                         zinc_chlnan,                         &
1470            &                         ln_slchlpicinc,                      &
1471            &                         zinc_chlpic,                         &
1472            &                         ln_slchldininc,                      &
1473            &                         zinc_chldin,                         &
1474            &                         zincper,                             &
1475            &                         rn_maxchlinc, ln_phytobal, zmld,     &
1476            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
1477            &                         phyt_avg_bkg, mld_max_bkg,           &
1478            &                         totalk_bkg,                          &
1479            &                         tracer_bkg, phyto2d_balinc )
1480#else
1481         CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
1482            &           'but not defined a biogeochemical model' )
1483#endif
1484
1485      ENDIF
1486
1487      IF ( ll_asmiau ) THEN
1488
1489         !--------------------------------------------------------------------
1490         ! Incremental Analysis Updating
1491         !--------------------------------------------------------------------
1492
1493         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1494
1495            it = kt - nit000 + 1
1496            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1497            ! note this is not a tendency so should not be divided by rdt
1498
1499            IF(lwp) THEN
1500               WRITE(numout,*) 
1501               WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
1502                  &  kt,' with IAU weight = ', pwgtiau(it)
1503               WRITE(numout,*) '~~~~~~~~~~~~'
1504            ENDIF
1505
1506            ! Update the biogeochemical variables
1507            ! Add directly to trn and trb, rather than to tra, because tra gets
1508            ! reset to zero at the start of trc_stp, called after this routine
1509#if defined key_medusa
1510            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1511               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1512               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1513                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1514               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1515                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1516            END WHERE
1517#elif defined key_hadocc
1518            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1519               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1520               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1521                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1522               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1523                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1524            END WHERE
1525#elif defined key_fabm
1526            WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1527               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
1528               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1529                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1530               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
1531                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1532            END WHERE
1533#endif
1534
1535            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1536            ! which is called at end of model run
1537         ENDIF
1538
1539      ELSEIF ( ll_asmdin ) THEN 
1540
1541         !--------------------------------------------------------------------
1542         ! Direct Initialization
1543         !--------------------------------------------------------------------
1544         
1545         IF ( kt == nitdin_r ) THEN
1546
1547            neuler = 0                    ! Force Euler forward step
1548
1549            ! Initialize the now fields with the background + increment
1550            ! Background currently is what the model is initialised with
1551            CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
1552               &           ' Background state is taken from model rather than background file' )
1553#if defined key_medusa
1554            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1555               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1556               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1557                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
1558               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1559            END WHERE
1560#elif defined key_hadocc
1561            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1562               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1563               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1564                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1)
1565               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1566            END WHERE
1567#elif defined key_fabm
1568            WHERE( phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1569               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
1570               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1571                  &                           phyto2d_balinc(:,:,:,jp_fabm0:jp_fabm1)
1572               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
1573            END WHERE
1574#endif
1575 
1576            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1577            ! which is called at end of model run
1578         ENDIF
1579         !
1580      ENDIF
1581      !
1582   END SUBROUTINE phyto2d_asm_inc
1583
1584   !!===========================================================================
1585   !!===========================================================================
1586   !!===========================================================================
1587
1588   SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1589      !!------------------------------------------------------------------------
1590      !!                    ***  ROUTINE phyto3d_asm_inc  ***
1591      !!         
1592      !! ** Purpose : Apply the profile chlorophyll assimilation increments.
1593      !!
1594      !! ** Method  : Calculate increments to state variables.
1595      !!              Direct initialization or Incremental Analysis Updating.
1596      !!
1597      !! ** Action  :
1598      !!------------------------------------------------------------------------
1599      INTEGER,  INTENT(IN) :: kt        ! Current time step
1600      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1601      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1602      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1603      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1604      !
1605      INTEGER                          :: ji, jj, jk    ! Loop counters
1606      INTEGER                          :: it            ! Index
1607      REAL(wp)                         :: zincwgt       ! IAU weight for timestep
1608      REAL(wp)                         :: zfrac_chn     ! Fraction of jpchn
1609      REAL(wp)                         :: zfrac_chd     ! Fraction of jpchd
1610      REAL(wp)                         :: zfrac_chl1    ! Fraction of jp_fabm_chl1
1611      REAL(wp)                         :: zfrac_chl2    ! Fraction of jp_fabm_chl2
1612      REAL(wp)                         :: zfrac_chl3    ! Fraction of jp_fabm_chl3
1613      REAL(wp)                         :: zfrac_chl4    ! Fraction of jp_fabm_chl4
1614      REAL(wp)                         :: zrat_phn_chn  ! jpphn:jpchn ratio
1615      REAL(wp)                         :: zrat_phd_chd  ! jpphd:jpchd ratio
1616      REAL(wp)                         :: zrat_pds_chd  ! jppds:jpchd ratio
1617      REAL(wp)                         :: zrat_p1c_chl1 ! jp_fabm_p1c:jp_fabm_chl1 ratio
1618      REAL(wp)                         :: zrat_p1n_chl1 ! jp_fabm_p1n:jp_fabm_chl1 ratio
1619      REAL(wp)                         :: zrat_p1p_chl1 ! jp_fabm_p1p:jp_fabm_chl1 ratio
1620      REAL(wp)                         :: zrat_p1s_chl1 ! jp_fabm_p1s:jp_fabm_chl1 ratio
1621      REAL(wp)                         :: zrat_p2c_chl2 ! jp_fabm_p2c:jp_fabm_chl2 ratio
1622      REAL(wp)                         :: zrat_p2n_chl2 ! jp_fabm_p2n:jp_fabm_chl2 ratio
1623      REAL(wp)                         :: zrat_p2p_chl2 ! jp_fabm_p2p:jp_fabm_chl2 ratio
1624      REAL(wp)                         :: zrat_p3c_chl3 ! jp_fabm_p3c:jp_fabm_chl3 ratio
1625      REAL(wp)                         :: zrat_p3n_chl3 ! jp_fabm_p3n:jp_fabm_chl3 ratio
1626      REAL(wp)                         :: zrat_p3p_chl3 ! jp_fabm_p3p:jp_fabm_chl3 ratio
1627      REAL(wp)                         :: zrat_p4c_chl4 ! jp_fabm_p4c:jp_fabm_chl4 ratio
1628      REAL(wp)                         :: zrat_p4n_chl4 ! jp_fabm_p4n:jp_fabm_chl4 ratio
1629      REAL(wp)                         :: zrat_p4p_chl4 ! jp_fabm_p4p:jp_fabm_chl4 ratio
1630      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc       ! Chlorophyll increments
1631      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl       ! Chlorophyll background
1632      !!------------------------------------------------------------------------
1633
1634      IF ( kt <= nit000 ) THEN
1635
1636         IF ( ln_plchltotinc ) THEN
1637            ! Convert log10(chlorophyll) increment back to a chlorophyll increment
1638            ! In order to transform logchl incs to chl incs, need to account for model
1639            ! background, cannot simply do 10^logchl_bkginc. Need to:
1640            ! 1) Add logchl inc to log10(background) to get log10(analysis)
1641            ! 2) Take 10^log10(analysis) to get analysis
1642            ! 3) Subtract background from analysis to get chl incs
1643            ! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
1644#if defined key_medusa
1645            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
1646#elif defined key_hadocc
1647            bkg_chl(:,:,:) = chl_bkg(:,:,:)
1648#elif defined key_fabm
1649            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl1) + &
1650               &             tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl2) + &
1651               &             tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl3) + &
1652               &             tracer_bkg(:,:,:,jp_fabm_m1+jp_fabm_chl4)
1653#endif
1654            DO jk = 1, jpk
1655               DO jj = 1, jpj
1656                  DO ji = 1, jpi
1657                     IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
1658                        chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
1659                        IF ( rn_maxchlinc > 0.0 ) THEN
1660                           chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
1661                        ENDIF
1662                     ELSE
1663                        chl_inc(ji,jj,jk) = 0.0
1664                     ENDIF
1665                  END DO
1666               END DO
1667            END DO
1668         ELSE IF ( ln_pchltotinc ) THEN
1669            DO jk = 1, jpk
1670               DO jj = 1, jpj
1671                  DO ji = 1, jpi
1672                     IF ( rn_maxchlinc > 0.0 ) THEN
1673                        chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
1674                     ELSE
1675                        chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
1676                     ENDIF
1677                  END DO
1678               END DO
1679            END DO
1680         ENDIF
1681
1682#if defined key_medusa
1683         ! Loop over each grid point partioning the increments based on existing ratios
1684         DO jk = 1, jpk
1685            DO jj = 1, jpj
1686               DO ji = 1, jpi
1687                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
1688                     zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
1689                     zfrac_chd = 1.0 - zfrac_chn
1690                     phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
1691                     phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
1692                     zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
1693                     zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
1694                     zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
1695                     phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
1696                     phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
1697                     phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
1698                  ENDIF
1699               END DO
1700            END DO
1701         END DO
1702#elif defined key_hadocc
1703         phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
1704#elif defined key_fabm
1705         ! Loop over each grid point partioning the increments based on existing ratios
1706         DO jk = 1, jpk
1707            DO jj = 1, jpj
1708               DO ji = 1, jpi
1709                  IF ( ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) > 0.0 ) .AND. &
1710                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) > 0.0 ) .AND. &
1711                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) > 0.0 ) .AND. &
1712                     & ( tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) > 0.0 ) ) THEN
1713                     zfrac_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) / bkg_chl(ji,jj,jk)
1714                     zfrac_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) / bkg_chl(ji,jj,jk)
1715                     zfrac_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) / bkg_chl(ji,jj,jk)
1716                     zfrac_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) / bkg_chl(ji,jj,jk)
1717                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) = chl_inc(ji,jj,jk) * zfrac_chl1
1718                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) = chl_inc(ji,jj,jk) * zfrac_chl2
1719                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) = chl_inc(ji,jj,jk) * zfrac_chl3
1720                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) = chl_inc(ji,jj,jk) * zfrac_chl4
1721                     zrat_p1c_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
1722                     zrat_p1n_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
1723                     zrat_p1p_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
1724                     zrat_p1s_chl1 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1)
1725                     zrat_p2c_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2)
1726                     zrat_p2n_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2)
1727                     zrat_p2p_chl2 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2)
1728                     zrat_p3c_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3)
1729                     zrat_p3n_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3)
1730                     zrat_p3p_chl3 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3)
1731                     zrat_p4c_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4)
1732                     zrat_p4n_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4)
1733                     zrat_p4p_chl4 = tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) / tracer_bkg(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4)
1734                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1c) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) * zrat_p1c_chl1
1735                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1n) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) * zrat_p1n_chl1
1736                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1p) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) * zrat_p1p_chl1
1737                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p1s) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl1) * zrat_p1s_chl1
1738                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2c) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) * zrat_p2c_chl2
1739                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2n) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) * zrat_p2n_chl2
1740                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p2p) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl2) * zrat_p2p_chl2
1741                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3c) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) * zrat_p3c_chl3
1742                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3n) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) * zrat_p3n_chl3
1743                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p3p) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl3) * zrat_p3p_chl3
1744                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4c) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) * zrat_p4c_chl4
1745                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4n) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) * zrat_p4n_chl4
1746                     phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_p4p) = phyto3d_balinc(ji,jj,jk,jp_fabm_m1+jp_fabm_chl4) * zrat_p4p_chl4
1747                  ENDIF
1748               END DO
1749            END DO
1750         END DO
1751#else
1752         CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
1753            &           'but not defined a biogeochemical model' )
1754#endif
1755
1756      ENDIF
1757
1758      IF ( ll_asmiau ) THEN
1759
1760         !--------------------------------------------------------------------
1761         ! Incremental Analysis Updating
1762         !--------------------------------------------------------------------
1763
1764         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1765
1766            it = kt - nit000 + 1
1767            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1768            ! note this is not a tendency so should not be divided by rdt
1769
1770            IF(lwp) THEN
1771               WRITE(numout,*) 
1772               WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
1773                  &  kt,' with IAU weight = ', pwgtiau(it)
1774               WRITE(numout,*) '~~~~~~~~~~~~'
1775            ENDIF
1776
1777            ! Update the biogeochemical variables
1778            ! Add directly to trn and trb, rather than to tra, because tra gets
1779            ! reset to zero at the start of trc_stp, called after this routine
1780#if defined key_medusa
1781            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1782               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1783               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1784                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1785               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1786                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1787            END WHERE
1788#elif defined key_hadocc
1789            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1790               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1791               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1792                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1793               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1794                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1795            END WHERE
1796#elif defined key_fabm
1797            WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1798               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
1799               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1800                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1801               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
1802                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
1803            END WHERE
1804#endif
1805
1806            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1807            ! which is called at end of model run
1808         ENDIF
1809
1810      ELSEIF ( ll_asmdin ) THEN 
1811
1812         !--------------------------------------------------------------------
1813         ! Direct Initialization
1814         !--------------------------------------------------------------------
1815         
1816         IF ( kt == nitdin_r ) THEN
1817
1818            neuler = 0                    ! Force Euler forward step
1819
1820            ! Initialize the now fields with the background + increment
1821            ! Background currently is what the model is initialised with
1822            CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
1823               &           ' Background state is taken from model rather than background file' )
1824#if defined key_medusa
1825            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1826               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1827               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1828                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
1829               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1830            END WHERE
1831#elif defined key_hadocc
1832            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1833               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1834               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1835                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1)
1836               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1837            END WHERE
1838#elif defined key_fabm
1839            WHERE( phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
1840               &   trn(:,:,:,jp_fabm0:jp_fabm1) + phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
1841               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
1842                  &                           phyto3d_balinc(:,:,:,jp_fabm0:jp_fabm1)
1843               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
1844            END WHERE
1845#endif
1846 
1847            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1848            ! which is called at end of model run
1849         ENDIF
1850         !
1851      ENDIF
1852      !
1853   END SUBROUTINE phyto3d_asm_inc
1854
1855   !!===========================================================================
1856   !!===========================================================================
1857   !!===========================================================================
1858
1859   SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1860      &                     ll_trainc, pt_bkginc, ps_bkginc )
1861      !!------------------------------------------------------------------------
1862      !!                    ***  ROUTINE pco2_asm_inc  ***
1863      !!         
1864      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
1865      !!
1866      !! ** Method  : Calculate increments to state variables using carbon
1867      !!              balancing.
1868      !!              Direct initialization or Incremental Analysis Updating.
1869      !!
1870      !! ** Action  :
1871      !!------------------------------------------------------------------------
1872      INTEGER, INTENT(IN)                   :: kt           ! Current time step
1873      LOGICAL, INTENT(IN)                   :: ll_asmdin    ! Flag for direct initialisation
1874      LOGICAL, INTENT(IN)                   :: ll_asmiau    ! Flag for incremental analysis update
1875      INTEGER, INTENT(IN)                   :: kcycper      ! Dimension of pwgtiau
1876      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1877      LOGICAL, INTENT(IN)                   :: ll_trainc    ! Flag for T/S increments
1878      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1879      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1880      !
1881      INTEGER                               :: ji, jj, jk   ! Loop counters
1882      INTEGER                               :: it           ! Index
1883      INTEGER                               :: jkmax        ! Index of mixed layer depth
1884      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
1885      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
1886      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
1887      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
1888      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
1889      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
1890      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
1891      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
1892      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
1893
1894      ! Coefficients for fCO2 to pCO2 conversion
1895      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
1896      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
1897      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
1898      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
1899      REAL(wp)                              :: zcoef_fco2_5 = 2.0
1900      REAL(wp)                              :: zcoef_fco2_6 = 57.7
1901      REAL(wp)                              :: zcoef_fco2_7 = 0.118
1902      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
1903      !!------------------------------------------------------------------------
1904
1905      IF ( kt <= nit000 ) THEN
1906
1907         !----------------------------------------------------------------------
1908         ! DIC and alkalinity balancing
1909         !----------------------------------------------------------------------
1910
1911         IF ( ln_sfco2inc ) THEN
1912#if defined key_medusa && defined key_roam
1913            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
1914            patm(1) = 1.0
1915            phyd(1) = 0.0
1916            DO jj = 1, jpj
1917               DO ji = 1, jpi
1918                  CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
1919               END DO
1920            END DO
1921#else
1922            ! If assimilating fCO2, then convert to pCO2 using temperature
1923            ! See flux_gas.F90 within HadOCC for details of calculation
1924            pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) /                                                            &
1925               &                    EXP((zcoef_fco2_1                                                            + &
1926               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
1927               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
1928               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
1929               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
1930               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
1931#endif
1932         ELSE
1933            pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
1934         ENDIF
1935
1936         ! Account for physics assimilation if required
1937         IF ( ll_trainc ) THEN
1938            tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
1939            sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
1940         ELSE
1941            tem_bkg_temp(:,:) = tsn(:,:,1,1)
1942            sal_bkg_temp(:,:) = tsn(:,:,1,2)
1943         ENDIF
1944
1945#if defined key_medusa
1946         ! Account for phytoplankton balancing if required
1947         IF ( ln_phytobal ) THEN
1948            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
1949            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
1950         ELSE
1951            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1952            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1953         ENDIF
1954
1955         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1956            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1957            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1958
1959#elif defined key_hadocc
1960         ! Account for phytoplankton balancing if required
1961         IF ( ln_phytobal ) THEN
1962            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
1963            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
1964         ELSE
1965            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1966            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1967         ENDIF
1968
1969         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1970            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1971            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1972
1973#elif defined key_fabm
1974         ! Account for phytoplankton balancing if required
1975         IF ( ln_phytobal ) THEN
1976            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_o3c) + phyto2d_balinc(:,:,1,jp_fabm_m1+jp_fabm_o3c)
1977            alk_bkg_temp(:,:) = totalk_bkg(:,:,1)             + phyto2d_balinc(:,:,1,jp_fabm_m1+jp_fabm_o3ba)
1978         ELSE
1979            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_fabm_m1+jp_fabm_o3c)
1980            alk_bkg_temp(:,:) = totalk_bkg(:,:,1)
1981         ENDIF
1982
1983         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1984            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1985            &               pco2_balinc(:,:,1,jp_fabm_m1+jp_fabm_o3c), pco2_balinc(:,:,1,jp_fabm_m1+jp_fabm_o3ba) )
1986
1987#else
1988         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1989            &           'but not defined a biogeochemical model' )
1990#endif
1991
1992         ! Select mixed layer
1993         IF ( ll_asmdin ) THEN
1994#if defined key_hadocc || defined key_medusa || defined key_fabm
1995            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1996               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1997            zmld(:,:) = mld_max_bkg(:,:)
1998#else
1999            CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
2000#endif
2001         ELSE
2002            SELECT CASE( mld_choice_bgc )
2003            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
2004               zmld(:,:) = hmld(:,:)
2005            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
2006               zmld(:,:) = hmlp(:,:)
2007            CASE ( 3 )                   ! Kara MLD [Interpolated]
2008#if defined key_karaml
2009               IF ( ln_kara ) THEN
2010                  zmld(:,:) = hmld_kara(:,:)
2011               ELSE
2012                  CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
2013                     &           ' but ln_kara=.false.' )
2014               ENDIF
2015#else
2016               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
2017                  &           ' but is not defined' )
2018#endif
2019            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
2020               zmld(:,:) = hmld_tref(:,:)
2021            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
2022               zmld(:,:) = hmlpt(:,:)
2023            END SELECT
2024         ENDIF
2025
2026         ! Propagate through mixed layer
2027         DO jj = 1, jpj
2028            DO ji = 1, jpi
2029               !
2030               jkmax = jpk-1
2031               DO jk = jpk-1, 1, -1
2032                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
2033                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
2034                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
2035                     jkmax = jk
2036                  ENDIF
2037               END DO
2038               !
2039#if defined key_top
2040               DO jk = 2, jkmax
2041                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
2042               END DO
2043#endif
2044               !
2045            END DO
2046         END DO
2047
2048      ENDIF
2049
2050      IF ( ll_asmiau ) THEN
2051
2052         !--------------------------------------------------------------------
2053         ! Incremental Analysis Updating
2054         !--------------------------------------------------------------------
2055
2056         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
2057
2058            it = kt - nit000 + 1
2059            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
2060            ! note this is not a tendency so should not be divided by rdt
2061
2062            IF(lwp) THEN
2063               WRITE(numout,*) 
2064               IF ( ln_spco2inc ) THEN
2065                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
2066                     &  kt,' with IAU weight = ', pwgtiau(it)
2067               ELSE IF ( ln_sfco2inc ) THEN
2068                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
2069                     &  kt,' with IAU weight = ', pwgtiau(it)
2070               ENDIF
2071               WRITE(numout,*) '~~~~~~~~~~~~'
2072            ENDIF
2073
2074            ! Update the biogeochemical variables
2075            ! Add directly to trn and trb, rather than to tra, because tra gets
2076            ! reset to zero at the start of trc_stp, called after this routine
2077#if defined key_medusa
2078            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2079               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
2080               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2081                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2082               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
2083                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2084            END WHERE
2085#elif defined key_hadocc
2086            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
2087               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
2088               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
2089                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
2090               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
2091                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
2092            END WHERE
2093#elif defined key_fabm
2094            WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
2095               &   trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt > 0.0_wp )
2096               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
2097                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
2098               trb(:,:,:,jp_fabm0:jp_fabm1) = trb(:,:,:,jp_fabm0:jp_fabm1) + &
2099                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) * zincwgt
2100            END WHERE
2101#endif
2102
2103            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
2104            ! which is called at end of model run
2105
2106         ENDIF
2107
2108      ELSEIF ( ll_asmdin ) THEN 
2109
2110         !--------------------------------------------------------------------
2111         ! Direct Initialization
2112         !--------------------------------------------------------------------
2113         
2114         IF ( kt == nitdin_r ) THEN
2115
2116            neuler = 0                    ! Force Euler forward step
2117
2118            ! Initialize the now fields with the background + increment
2119            ! Background currently is what the model is initialised with
2120            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
2121               &           ' Background state is taken from model rather than background file' )
2122#if defined key_medusa
2123            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2124               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
2125               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2126                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1)
2127               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
2128            END WHERE
2129#elif defined key_hadocc
2130            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
2131               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
2132               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
2133                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
2134               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
2135            END WHERE
2136#elif defined key_fabm
2137            WHERE( pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp .OR. &
2138               &   trn(:,:,:,jp_fabm0:jp_fabm1) + pco2_balinc(:,:,:,jp_fabm0:jp_fabm1) > 0.0_wp )
2139               trn(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1) + &
2140                  &                           pco2_balinc(:,:,:,jp_fabm0:jp_fabm1)
2141               trb(:,:,:,jp_fabm0:jp_fabm1) = trn(:,:,:,jp_fabm0:jp_fabm1)
2142            END WHERE
2143#endif
2144 
2145            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
2146            ! which is called at end of model run
2147         ENDIF
2148         !
2149      ENDIF
2150      !
2151   END SUBROUTINE pco2_asm_inc
2152
2153   !!===========================================================================
2154   !!===========================================================================
2155   !!===========================================================================
2156
2157   SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
2158      &                   ll_trainc, pt_bkginc, ps_bkginc )
2159      !!------------------------------------------------------------------------
2160      !!                    ***  ROUTINE ph_asm_inc  ***
2161      !!         
2162      !! ** Purpose : Apply the pH assimilation increments.
2163      !!
2164      !! ** Method  : Calculate increments to DIC and alkalinity from pH
2165      !!              Use a similar approach to the pCO2 scheme
2166      !!              Direct initialization or Incremental Analysis Updating.
2167      !!
2168      !! ** Action  :
2169      !!------------------------------------------------------------------------
2170      INTEGER,                          INTENT(IN) :: kt        ! Current time step
2171      LOGICAL,                          INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
2172      LOGICAL,                          INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
2173      INTEGER,                          INTENT(IN) :: kcycper   ! Dimension of pwgtiau
2174      REAL(wp), DIMENSION(kcycper),     INTENT(IN) :: pwgtiau   ! IAU weights
2175      LOGICAL,                          INTENT(IN) :: ll_trainc ! Flag for T/S increments
2176      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
2177      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
2178     
2179      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
2180      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
2181      REAL(wp)                         :: PH_OUT         ! pH from pH calculation
2182      REAL(wp)                         :: ph_bkg         ! pH from pH calculation
2183      REAL(wp)                         :: ph_dic         ! pH from pH calculation
2184      REAL(wp)                         :: ph_alk         ! pH from pH calculation
2185      REAL(wp)                         :: dph_ddic       ! Change in pH wrt DIC
2186      REAL(wp)                         :: dph_dalk       ! Change in pH wrt alk
2187      REAL(wp)                         :: weight         ! Increment weighting
2188      REAL(wp)                         :: zincwgt        ! IAU weight for current time step
2189      REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp   ! DIC background
2190      REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp   ! Alkalinity background
2191      REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp   ! DIN background
2192      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp   ! Silicate background
2193      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp   ! Temperature background
2194      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp   ! Salinity background
2195      REAL(wp), DIMENSION(19)          :: dummy_out      ! Dummy variables
2196      INTEGER                          :: ji, jj, jk, jx ! Loop counters
2197      INTEGER                          :: it             ! Index
2198      !!------------------------------------------------------------------------
2199
2200#if ! defined key_medusa
2201      CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
2202#else
2203
2204      IF ( kt <= nit000 ) THEN
2205
2206         !----------------------------------------------------------------------
2207         ! DIC and alkalinity balancing
2208         !----------------------------------------------------------------------
2209
2210         ! Account for physics assimilation if required
2211         IF ( ll_trainc ) THEN
2212            tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
2213            sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
2214         ELSE
2215            tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
2216            sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
2217         ENDIF
2218
2219         ! Account for phytoplankton balancing if required
2220         IF ( ln_phytobal ) THEN
2221            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
2222            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
2223            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
2224            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
2225         ELSE
2226            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
2227            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
2228            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
2229            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
2230         ENDIF
2231         
2232         ! Account for pCO2 balancing if required
2233         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2234            dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
2235            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
2236         ENDIF
2237         
2238         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
2239         ! This requires three calls to the MOCSY carbonate package
2240         ! One to find pH at background DIC and alkalinity
2241         ! One to find pH when DIC is increased by 10
2242         ! One to find pH when alkalinity is increased by 10
2243         ! Then calculate the gradients
2244
2245         DO jk = 1, jpk
2246            DO jj = 1, jpj
2247               DO ji = 1, jpi
2248
2249                  IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
2250
2251                     DO jx = 1, 3
2252                        IF ( jx == 1 ) THEN
2253                           DIC_IN = dic_bkg_temp(ji,jj,jk)
2254                           ALK_IN = alk_bkg_temp(ji,jj,jk)
2255                        ELSE IF ( jx == 2 ) THEN
2256                           DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
2257                           ALK_IN = alk_bkg_temp(ji,jj,jk)
2258                        ELSE IF ( jx == 3 ) THEN
2259                           DIC_IN = dic_bkg_temp(ji,jj,jk)
2260                           ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
2261                        ENDIF
2262
2263                        CALL mocsy_interface(tsn(ji,jj,jk,jp_tem),         & ! IN:  temperature
2264                           &                 tsn(ji,jj,jk,jp_sal),         & !      salinity
2265                           &                 ALK_IN,                       & !      alkalinity
2266                           &                 DIC_IN,                       & !      DIC
2267                           &                 sil_bkg_temp(ji,jj,jk),       & !      silicate
2268                           &                 din_bkg_temp(ji,jj,jk)/16.0,  & !      phosphate (Redfield from DIN)
2269                           &                 1.0,                          & !      atmospheric pressure (dummy)
2270                           &                 fsdept(ji,jj,jk),             & !      depth
2271                           &                 gphit(ji,jj),                 & !      latitude
2272                           &                 1.0,                          & !      gas transfer (dummy)
2273                           &                 1.0,                          & !      atmospheric xCO2 (dummy)
2274                           &                 1,                            & !      number of points
2275                           &                 PH_OUT,                       & ! OUT: pH
2276                           &                 dummy_out(1),  dummy_out(2),  & !      stuff we don't care about
2277                           &                 dummy_out(3),  dummy_out(4),  &
2278                           &                 dummy_out(5),  dummy_out(6),  &
2279                           &                 dummy_out(7),  dummy_out(8),  &
2280                           &                 dummy_out(9),  dummy_out(10), &
2281                           &                 dummy_out(11), dummy_out(12), &
2282                           &                 dummy_out(13), dummy_out(14), &
2283                           &                 dummy_out(15), dummy_out(16), &
2284                           &                 dummy_out(17), dummy_out(18), &
2285                           &                 dummy_out(19) )
2286
2287                        IF ( jx == 1 ) THEN
2288                           ph_bkg = PH_OUT
2289                        ELSE IF ( jx == 2 ) THEN
2290                           ph_dic = PH_OUT
2291                        ELSE IF ( jx == 3 ) THEN
2292                           ph_alk = PH_OUT
2293                        ENDIF
2294                     END DO
2295
2296                     dph_ddic = (ph_dic - ph_bkg) / zsearch
2297                     dph_dalk = (ph_alk - ph_bkg) / zsearch
2298                     weight   = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) )
2299
2300                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
2301                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
2302                     
2303                  ENDIF
2304                 
2305               END DO
2306            END DO
2307         END DO
2308
2309      ENDIF
2310     
2311      IF ( ll_asmiau ) THEN
2312
2313         !--------------------------------------------------------------------
2314         ! Incremental Analysis Updating
2315         !--------------------------------------------------------------------
2316
2317         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
2318
2319            it = kt - nit000 + 1
2320            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
2321            ! note this is not a tendency so should not be divided by rdt
2322
2323            IF(lwp) THEN
2324               WRITE(numout,*) 
2325               WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
2326                  &  kt,' with IAU weight = ', pwgtiau(it)
2327               WRITE(numout,*) '~~~~~~~~~~~~'
2328            ENDIF
2329
2330            ! Update the biogeochemical variables
2331            ! Add directly to trn and trb, rather than to tra, because tra gets
2332            ! reset to zero at the start of trc_stp, called after this routine
2333            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2334               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
2335               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2336                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2337               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
2338                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
2339            END WHERE
2340
2341            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
2342            ! which is called at end of model run
2343
2344         ENDIF
2345
2346      ELSEIF ( ll_asmdin ) THEN 
2347
2348         !--------------------------------------------------------------------
2349         ! Direct Initialization
2350         !--------------------------------------------------------------------
2351         
2352         IF ( kt == nitdin_r ) THEN
2353
2354            neuler = 0                    ! Force Euler forward step
2355
2356            ! Initialize the now fields with the background + increment
2357            ! Background currently is what the model is initialised with
2358            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', &
2359               &           ' Background state is taken from model rather than background file' )
2360            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
2361               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
2362               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
2363                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1)
2364               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
2365            END WHERE
2366 
2367            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
2368            ! which is called at end of model run
2369         ENDIF
2370         !
2371      ENDIF
2372#endif     
2373      !
2374   END SUBROUTINE ph_asm_inc
2375
2376   !!===========================================================================
2377   !!===========================================================================
2378   !!===========================================================================
2379
2380   SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
2381      !!----------------------------------------------------------------------
2382      !!                    ***  ROUTINE dyn_asm_inc  ***
2383      !!         
2384      !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
2385      !!
2386      !! ** Method  : Direct initialization or Incremental Analysis Updating.
2387      !!
2388      !! ** Action  :
2389      !!----------------------------------------------------------------------
2390      INTEGER,  INTENT(IN) :: kt        ! Current time step
2391      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
2392      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
2393      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
2394      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
2395      !
2396      INTEGER  :: it              ! Index
2397      REAL(wp) :: zincwgt         ! IAU weight for current time step
2398      !!----------------------------------------------------------------------
2399
2400      IF ( kt <= nit000 ) THEN
2401
2402         !----------------------------------------------------------------------
2403         ! Remove any other balancing increments
2404         !----------------------------------------------------------------------
2405
2406         IF ( ln_pno3inc ) THEN
2407#if defined key_hadocc || defined key_medusa || defined key_fabm
2408#if defined key_hadocc
2409            it = jp_had_nut
2410#elif defined key_medusa
2411            it = jpdin
2412#elif defined key_fabm
2413            it = jp_fabm_n3n
2414#endif
2415            IF ( ln_phytobal ) THEN
2416               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2417            ENDIF
2418            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2419               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2420            ENDIF
2421            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2422               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2423            ENDIF
2424            IF ( ln_pphinc ) THEN
2425               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2426            ENDIF
2427#else
2428            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2429#endif
2430         ENDIF
2431
2432         IF ( ln_psi4inc ) THEN
2433#if defined key_medusa || defined key_fabm
2434#if defined key_medusa
2435            it = jpsil
2436#elif defined key_fabm
2437            it = jp_fabm_n5s
2438#endif
2439            IF ( ln_phytobal ) THEN
2440               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2441            ENDIF
2442            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2443               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2444            ENDIF
2445            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2446               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2447            ENDIF
2448            IF ( ln_pphinc ) THEN
2449               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2450            ENDIF
2451#else
2452            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2453#endif
2454         ENDIF
2455
2456         IF ( ln_ppo4inc ) THEN
2457#if defined key_fabm
2458            it = jp_fabm_n1p
2459
2460            IF ( ln_phytobal ) THEN
2461               ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2462            ENDIF
2463            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2464               ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2465            ENDIF
2466            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2467               ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2468            ENDIF
2469            IF ( ln_pphinc ) THEN
2470               ppo4_bkginc(:,:,:) = ppo4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2471            ENDIF
2472#else
2473            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2474#endif
2475         ENDIF
2476
2477         IF ( ln_pdicinc ) THEN
2478#if defined key_hadocc || defined key_medusa || defined key_fabm
2479#if defined key_hadocc
2480            it = jp_had_dic
2481#elif defined key_medusa
2482            it = jpdic
2483#elif defined key_fabm
2484            it = jp_fabm_o3c
2485#endif
2486            IF ( ln_phytobal ) THEN
2487               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2488            ENDIF
2489            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2490               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2491            ENDIF
2492            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2493               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2494            ENDIF
2495            IF ( ln_pphinc ) THEN
2496               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2497            ENDIF
2498#else
2499            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2500#endif
2501         ENDIF
2502
2503         IF ( ln_palkinc ) THEN
2504#if defined key_hadocc || defined key_medusa || defined key_fabm
2505#if defined key_hadocc
2506            it = jp_had_alk
2507#elif defined key_medusa
2508            it = jpalk
2509#elif defined key_fabm
2510            it = jp_fabm_o3ba
2511#endif
2512            IF ( ln_phytobal ) THEN
2513               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2514            ENDIF
2515            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2516               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2517            ENDIF
2518            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2519               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2520            ENDIF
2521            IF ( ln_pphinc ) THEN
2522               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2523            ENDIF
2524#else
2525            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2526#endif
2527         ENDIF
2528
2529         IF ( ln_po2inc ) THEN
2530#if defined key_medusa || defined key_fabm
2531#if defined key_medusa
2532            it = jpoxy
2533#elif defined key_fabm
2534            it = jp_fabm_o2o
2535#endif
2536            IF ( ln_phytobal ) THEN
2537               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
2538            ENDIF
2539            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
2540               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
2541            ENDIF
2542            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
2543               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
2544            ENDIF
2545            IF ( ln_pphinc ) THEN
2546               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
2547            ENDIF
2548#else
2549            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2550#endif
2551         ENDIF
2552
2553      ENDIF
2554
2555      IF ( ll_asmiau ) THEN
2556
2557         !--------------------------------------------------------------------
2558         ! Incremental Analysis Updating
2559         !--------------------------------------------------------------------
2560
2561         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
2562
2563            it = kt - nit000 + 1
2564            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
2565            ! note this is not a tendency so should not be divided by rdt
2566
2567            IF(lwp) THEN
2568               WRITE(numout,*) 
2569               WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
2570                  &  kt,' with IAU weight = ', pwgtiau(it)
2571               WRITE(numout,*) '~~~~~~~~~~~~'
2572            ENDIF
2573
2574            ! Update the 3D BGC variables
2575            ! Add directly to trn and trb, rather than to tra, because tra gets
2576            ! reset to zero at the start of trc_stp, called after this routine
2577            ! Don't apply increments if they'll take concentrations negative
2578
2579            IF ( ln_pno3inc ) THEN
2580#if defined key_hadocc
2581               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2582                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2583                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2584                  trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2585               END WHERE
2586#elif defined key_medusa
2587               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2588                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2589                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2590                  trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2591               END WHERE
2592#elif defined key_fabm
2593               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2594                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2595                  trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
2596                  trb(:,:,:,jp_fabm_m1+jp_fabm_n3n) = trb(:,:,:,jp_fabm_m1+jp_fabm_n3n) + pno3_bkginc(:,:,:) * zincwgt
2597               END WHERE
2598#else
2599               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2600#endif
2601            ENDIF
2602
2603            IF ( ln_psi4inc ) THEN
2604#if defined key_medusa
2605               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2606                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2607                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2608                  trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2609               END WHERE
2610#elif defined key_fabm
2611               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2612                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2613                  trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
2614                  trb(:,:,:,jp_fabm_m1+jp_fabm_n5s) = trb(:,:,:,jp_fabm_m1+jp_fabm_n5s) + psi4_bkginc(:,:,:) * zincwgt
2615               END WHERE
2616#else
2617               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2618#endif
2619            ENDIF
2620
2621            IF ( ln_ppo4inc ) THEN
2622#if defined key_fabm
2623               WHERE( ppo4_bkginc(:,:,:) > 0.0_wp .OR. &
2624                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) + ppo4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2625                  trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) + ppo4_bkginc(:,:,:) * zincwgt
2626                  trb(:,:,:,jp_fabm_m1+jp_fabm_n1p) = trb(:,:,:,jp_fabm_m1+jp_fabm_n1p) + ppo4_bkginc(:,:,:) * zincwgt
2627               END WHERE
2628#else
2629               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2630#endif
2631            ENDIF
2632
2633            IF ( ln_pdicinc ) THEN
2634#if defined key_hadocc
2635               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2636                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2637                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2638                  trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2639               END WHERE
2640#elif defined key_medusa
2641               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2642                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2643                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2644                  trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2645               END WHERE
2646#elif defined key_fabm
2647               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2648                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2649                  trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
2650                  trb(:,:,:,jp_fabm_m1+jp_fabm_o3c) = trb(:,:,:,jp_fabm_m1+jp_fabm_o3c) + pdic_bkginc(:,:,:) * zincwgt
2651               END WHERE
2652#else
2653               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2654#endif
2655            ENDIF
2656
2657            IF ( ln_palkinc ) THEN
2658#if defined key_hadocc
2659               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2660                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2661                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2662                  trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2663               END WHERE
2664#elif defined key_medusa
2665               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2666                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2667                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2668                  trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2669               END WHERE
2670#elif defined key_fabm
2671               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2672                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2673                  trn(:,:,:,jp_fabm_m1+jp_fabm_o3ba) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
2674                  trb(:,:,:,jp_fabm_m1+jp_fabm_o3ba) = trb(:,:,:,jp_fabm_m1+jp_fabm_o3ba) + palk_bkginc(:,:,:) * zincwgt
2675               END WHERE
2676#else
2677               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2678#endif
2679            ENDIF
2680
2681            IF ( ln_po2inc ) THEN
2682#if defined key_medusa
2683               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2684                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
2685                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2686                  trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2687               END WHERE
2688#elif defined key_fabm
2689               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2690                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
2691                  trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
2692                  trb(:,:,:,jp_fabm_m1+jp_fabm_o2o) = trb(:,:,:,jp_fabm_m1+jp_fabm_o2o) + po2_bkginc(:,:,:) * zincwgt
2693               END WHERE
2694#else
2695               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2696#endif
2697            ENDIF
2698           
2699            IF ( kt == nitiaufin_r ) THEN
2700               IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2701               IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2702               IF ( ln_ppo4inc ) DEALLOCATE( ppo4_bkginc )
2703               IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2704               IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2705               IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
2706            ENDIF
2707
2708         ENDIF
2709
2710      ELSEIF ( ll_asmdin ) THEN 
2711
2712         !--------------------------------------------------------------------
2713         ! Direct Initialization
2714         !--------------------------------------------------------------------
2715         
2716         IF ( kt == nitdin_r ) THEN
2717
2718            neuler = 0                    ! Force Euler forward step
2719
2720            ! Initialize the now fields with the background + increment
2721            ! Background currently is what the model is initialised with
2722            CALL ctl_warn( ' Doing direct initialisation with 3D BGC assimilation', &
2723               &           ' Background state is taken from model rather than background file' )
2724
2725            IF ( ln_pno3inc ) THEN
2726#if defined key_hadocc
2727               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2728                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
2729                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
2730                  trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
2731               END WHERE
2732#elif defined key_medusa
2733               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2734                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
2735                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
2736                  trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
2737               END WHERE
2738#elif defined key_fabm
2739               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2740                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + pno3_bkginc(:,:,:) > 0.0_wp )
2741                  trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n) + pno3_bkginc(:,:,:)
2742                  trb(:,:,:,jp_fabm_m1+jp_fabm_n3n) = trn(:,:,:,jp_fabm_m1+jp_fabm_n3n)
2743               END WHERE
2744#else
2745               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2746#endif
2747            ENDIF
2748
2749            IF ( ln_psi4inc ) THEN
2750#if defined key_medusa
2751               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2752                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
2753                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
2754                  trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
2755               END WHERE
2756#elif defined key_fabm
2757               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2758                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) + psi4_bkginc(:,:,:) > 0.0_wp )
2759                  trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s) + psi4_bkginc(:,:,:)
2760                  trb(:,:,:,jp_fabm_m1+jp_fabm_n5s) = trn(:,:,:,jp_fabm_m1+jp_fabm_n5s)
2761               END WHERE
2762#else
2763               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2764#endif
2765            ENDIF
2766
2767            IF ( ln_ppo4inc ) THEN
2768#if defined key_fabm
2769               WHERE( ppo4_bkginc(:,:,:) > 0.0_wp .OR. &
2770                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) + ppo4_bkginc(:,:,:) > 0.0_wp )
2771                  trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p) + ppo4_bkginc(:,:,:)
2772                  trb(:,:,:,jp_fabm_m1+jp_fabm_n1p) = trn(:,:,:,jp_fabm_m1+jp_fabm_n1p)
2773               END WHERE
2774#else
2775               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2776#endif
2777            ENDIF
2778
2779            IF ( ln_pdicinc ) THEN
2780#if defined key_hadocc
2781               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2782                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
2783                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
2784                  trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
2785               END WHERE
2786#elif defined key_medusa
2787               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2788                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
2789                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
2790                  trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
2791               END WHERE
2792#elif defined key_fabm
2793               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2794                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) + pdic_bkginc(:,:,:) > 0.0_wp )
2795                  trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c) + pdic_bkginc(:,:,:)
2796                  trb(:,:,:,jp_fabm_m1+jp_fabm_o3c) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3c)
2797               END WHERE
2798#else
2799               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2800#endif
2801            ENDIF
2802
2803            IF ( ln_palkinc ) THEN
2804#if defined key_hadocc
2805               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2806                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
2807                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
2808                  trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
2809               END WHERE
2810#elif defined key_medusa
2811               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2812                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
2813                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
2814                  trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
2815               END WHERE
2816#elif defined key_fabm
2817               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2818                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_o3ba) + palk_bkginc(:,:,:) > 0.0_wp )
2819                  trn(:,:,:,jp_fabm_m1+jp_fabm_o3ba) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ba) + palk_bkginc(:,:,:)
2820                  trb(:,:,:,jp_fabm_m1+jp_fabm_o3ba) = trn(:,:,:,jp_fabm_m1+jp_fabm_o3ba)
2821               END WHERE
2822#else
2823               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2824#endif
2825            ENDIF
2826
2827            IF ( ln_po2inc ) THEN
2828#if defined key_medusa
2829               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2830                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
2831                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
2832                  trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
2833               END WHERE
2834#elif defined key_fabm
2835               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2836                  &   trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) + po2_bkginc(:,:,:) > 0.0_wp )
2837                  trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o) + po2_bkginc(:,:,:)
2838                  trb(:,:,:,jp_fabm_m1+jp_fabm_o2o) = trn(:,:,:,jp_fabm_m1+jp_fabm_o2o)
2839               END WHERE
2840#else
2841               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2842#endif
2843            ENDIF
2844 
2845            IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2846            IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2847            IF ( ln_ppo4inc ) DEALLOCATE( ppo4_bkginc )
2848            IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2849            IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2850            IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
2851         ENDIF
2852         !
2853      ENDIF
2854      !
2855   END SUBROUTINE bgc3d_asm_inc
2856
2857   !!===========================================================================
2858
2859END MODULE asmbgc
Note: See TracBrowser for help on using the repository browser.