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

Last change on this file since 13318 was 13318, checked in by dford, 5 months ago

Allow assimilation of PFT chlorophyll profiles. See Met Office utils ticket 346.

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