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

Last change on this file since 10624 was 10624, checked in by dford, 18 months ago

Couple of bug fixes.

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