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

Last change on this file since 10622 was 10622, checked in by dford, 19 months ago

Implement biogeochemistry assimilation for FABM-ERSEM.

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