source: branches/UKMO/AMM15_v3_6_STABLE_package_collate/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10728

Last change on this file since 10728 was 10728, checked in by dford, 22 months ago

Changes for assimilation of biogeochemical variables. See https://code.metoffice.gov.uk/trac/utils/ticket/174.

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