New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asmbgc.F90 in branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10149

Last change on this file since 10149 was 10149, checked in by frrh, 6 years ago

Met Office GMED ticket 379: Merged David Ford's MEDUSA assimilation changes
using command:

svn merge -r 10054:10141 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_v3

File size: 102.3 KB
Line 
1MODULE asmbgc
2   !!===========================================================================
3   !!                       ***  MODULE asmbgc  ***
4   !! Routines and declarations for biogeochemical assimilation
5   !!===========================================================================
6   !! History :  3.6  ! 2018-02 (D. Ford)  Adapted from existing modules
7   !!---------------------------------------------------------------------------
8   !! 'key_asminc'          : assimilation increment interface
9   !! 'key_top'             : passive tracers
10   !! 'key_hadocc'          : HadOCC model
11   !! 'key_medusa'          : MEDUSA model
12   !! 'key_roam'            : MEDUSA extras for carbonate cycle
13   !! 'key_karaml'          : Kara mixed layer depth
14   !!---------------------------------------------------------------------------
15   !! asm_bgc_check_options : check bgc assimilation options
16   !! asm_bgc_init_incs     : initialise bgc increments
17   !! asm_bgc_init_bkg      : initialise bgc background
18   !! asm_bgc_bal_wri       : write out bgc balancing increments
19   !! asm_bgc_bkg_wri       : write out bgc background
20   !! phyto2d_asm_inc       : apply the ocean colour increments
21   !! phyto3d_asm_inc       : apply the 3D phytoplankton increments
22   !! pco2_asm_inc          : apply the pCO2/fCO2 increments
23   !! ph_asm_inc            : apply the pH increments
24   !! bgc3d_asm_inc         : apply the generic 3D BGC increments
25   !!---------------------------------------------------------------------------
26   USE par_kind, ONLY:      & ! kind parameters
27      & wp
28   USE par_oce, ONLY:       & ! domain array sizes
29      & jpi,                &
30      & jpj,                &
31      & jpk
32   USE iom                    ! i/o
33   USE oce, ONLY:           & ! active tracer variables
34      & tsn
35   USE zdfmxl, ONLY :       & ! mixed layer depth
36#if defined key_karaml
37      & hmld_kara,          &
38      & ln_kara,            &
39#endif   
40      & hmld,               & 
41      & hmlp,               &
42      & hmlpt 
43   USE asmpar, ONLY:        & ! assimilation parameters
44      & c_asmbkg,           &
45      & c_asmbal,           &
46      & nitbkg_r,           &
47      & nitdin_r,           &
48      & nitiaustr_r,        &
49      & nitiaufin_r
50#if defined key_top
51   USE trc, ONLY:           & ! passive tracer variables
52      & trn,                &
53      & trb
54   USE par_trc, ONLY:       & ! passive tracer parameters
55      & jptra
56#endif
57#if defined key_medusa
58   USE asmphyto2dbal_medusa, ONLY: & ! phyto2d balancing for MEDUSA
59      & asm_phyto2d_bal_medusa
60   USE asmpco2bal, ONLY:    & ! pCO2 balancing for MEDUSA
61      & asm_pco2_bal
62   USE sms_medusa             ! MEDUSA parameters
63   USE par_medusa             ! MEDUSA parameters
64   USE mocsy_wrapper, ONLY: & ! MEDUSA carbonate system
65      & mocsy_interface
66   USE mocsy_mainmod, ONLY: & ! fCO2 to pCO2 conversion
67      & f2pCO2
68   USE sms_medusa, ONLY:    & ! MEDUSA diagnostic variables
69      & pgrow_avg,          &
70      & ploss_avg,          &
71      & phyt_avg,           &
72      & mld_max
73#elif defined key_hadocc
74   USE asmphyto2dbal_hadocc, ONLY: & ! phyto2d balancing for HadOCC
75      & asm_phyto2d_bal_hadocc
76   USE asmpco2bal, ONLY:    & ! pCO2 balancing for HadOCC
77      & asm_pco2_bal
78   USE par_hadocc             ! HadOCC parameters
79   USE had_bgc_const          ! HadOCC parameters
80   USE trc, ONLY:           & ! HadOCC diagnostic variables
81      & pgrow_avg,          &
82      & ploss_avg,          &
83      & phyt_avg,           &
84      & mld_max,            &
85      & HADOCC_CHL
86   USE had_bgc_const, ONLY: & ! HadOCC C:Chl ratio
87      & cchl_p
88#endif
89
90   IMPLICIT NONE
91   PRIVATE                   
92
93   PUBLIC  asm_bgc_check_options  ! called by asm_inc_init in asminc.F90
94   PUBLIC  asm_bgc_init_incs      ! called by asm_inc_init in asminc.F90
95   PRIVATE asm_bgc_read_incs_2d   ! called by asm_bgc_init_incs
96   PRIVATE asm_bgc_read_incs_3d   ! called by asm_bgc_init_incs
97   PUBLIC  asm_bgc_init_bkg       ! called by asm_inc_init in asminc.F90
98   PUBLIC  asm_bgc_bal_wri        ! called by nemo_gcm in nemogcm.F90
99   PUBLIC  asm_bgc_bkg_wri        ! called by asm_bkg_wri in asmbkg.F90
100   PUBLIC  phyto2d_asm_inc        ! called by bgc_asm_inc in asminc.F90
101   PUBLIC  phyto3d_asm_inc        ! called by bgc_asm_inc in asminc.F90
102   PUBLIC  pco2_asm_inc           ! called by bgc_asm_inc in asminc.F90
103   PUBLIC  ph_asm_inc             ! called by bgc_asm_inc in asminc.F90
104   PUBLIC  bgc3d_asm_inc          ! called by bgc_asm_inc in asminc.F90
105
106   LOGICAL, PUBLIC :: ln_balwri      = .FALSE. !: No output of balancing incs
107   LOGICAL, PUBLIC :: ln_phytobal    = .FALSE. !: No phytoplankton balancing
108   LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE. !: No surface total      log10(chlorophyll) increment
109   LOGICAL, PUBLIC :: ln_slchldiainc = .FALSE. !: No surface diatom     log10(chlorophyll) increment
110   LOGICAL, PUBLIC :: ln_slchlnoninc = .FALSE. !: No surface non-diatom log10(chlorophyll) increment
111   LOGICAL, PUBLIC :: ln_schltotinc  = .FALSE. !: No surface total      chlorophyll        increment
112   LOGICAL, PUBLIC :: ln_slphytotinc = .FALSE. !: No surface total      log10(phyto C)     increment
113   LOGICAL, PUBLIC :: ln_slphydiainc = .FALSE. !: No surface diatom     log10(phyto C)     increment
114   LOGICAL, PUBLIC :: ln_slphynoninc = .FALSE. !: No surface non-diatom log10(phyto C)     increment
115   LOGICAL, PUBLIC :: ln_spco2inc    = .FALSE. !: No surface pCO2                          increment
116   LOGICAL, PUBLIC :: ln_sfco2inc    = .FALSE. !: No surface fCO2                          increment
117   LOGICAL, PUBLIC :: ln_plchltotinc = .FALSE. !: No profile total      log10(chlorophyll) increment
118   LOGICAL, PUBLIC :: ln_pchltotinc  = .FALSE. !: No profile total      chlorophyll        increment
119   LOGICAL, PUBLIC :: ln_pno3inc     = .FALSE. !: No profile nitrate                       increment
120   LOGICAL, PUBLIC :: ln_psi4inc     = .FALSE. !: No profile silicate                      increment
121   LOGICAL, PUBLIC :: ln_pdicinc     = .FALSE. !: No profile dissolved inorganic carbon    increment
122   LOGICAL, PUBLIC :: ln_palkinc     = .FALSE. !: No profile alkalinity                    increment
123   LOGICAL, PUBLIC :: ln_pphinc      = .FALSE. !: No profile pH                            increment
124   LOGICAL, PUBLIC :: ln_po2inc      = .FALSE. !: No profile oxygen                        increment
125
126   INTEGER, PUBLIC :: mld_choice_bgc = 1 !: choice of mld for bgc assimilation
127                                         !  1) hmld      - Turbocline/mixing depth
128                                         !                 [W points]
129                                         !  2) hmlp      - Density criterion
130                                         !                 (0.01 kg/m^3 change from 10m)
131                                         !                 [W points]
132                                         !  3) hmld_kara - Kara MLD
133                                         !                 [Interpolated]
134                                         !  4) hmld_tref - Temperature criterion
135                                         !                 (0.2 K change from surface)
136                                         !                 [T points]
137                                         !  5) hmlpt     - Density criterion
138                                         !                 (0.01 kg/m^3 change from 10m)
139                                         !                 [T points]
140
141   REAL(wp), PUBLIC :: rn_maxchlinc = -999.0 !: maximum absolute non-log chlorophyll inc
142                                             ! <= 0 no maximum applied (switch turned off)
143                                             !  > 0 absolute chl inc capped at this value
144
145   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchltot_bkginc ! slchltot inc
146   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchldia_bkginc ! slchldia inc
147   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slchlnon_bkginc ! slchlnon inc
148   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: schltot_bkginc  ! schltot inc
149   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphytot_bkginc ! slphytot inc
150   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphydia_bkginc ! slphydia inc
151   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: slphynon_bkginc ! slphynon inc
152   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: sfco2_bkginc    ! sfco2 inc
153   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: spco2_bkginc    ! spco2 inc
154   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: plchltot_bkginc ! plchltot inc
155   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pchltot_bkginc  ! pchltot inc
156   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pno3_bkginc     ! pno3 inc
157   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: psi4_bkginc     ! psi4 inc
158   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pdic_bkginc     ! pdic inc
159   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: palk_bkginc     ! palk inc
160   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: pph_bkginc      ! pph inc
161   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: po2_bkginc      ! po2 inc
162   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto2d_balinc  ! Balancing incs from ocean colour
163   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: phyto3d_balinc  ! Balancing incs from p(l)chltot
164   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: pco2_balinc     ! Balancing incs from spco2/sfco2
165   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ph_balinc       ! Balancing incs from pph
166
167   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg   ! Background phyto growth
168   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg   ! Background phyto loss
169   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg    ! Background phyto conc
170   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg     ! Background max MLD
171   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg      ! Background tracer state
172   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: chl_bkg         ! Background chl
173   REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE :: cchl_p_bkg      ! Background c:chl
174
175#  include "domzgr_substitute.h90"
176
177CONTAINS
178
179   SUBROUTINE asm_bgc_check_options
180      !!------------------------------------------------------------------------
181      !!                    ***  ROUTINE asm_bgc_check_options  ***
182      !!
183      !! ** Purpose :   check validity of biogeochemical assimilation options
184      !!
185      !! ** Method  :   check settings of logicals
186      !!
187      !! ** Action  :   call ctl_stop/ctl_warn if required
188      !!
189      !! References :   asm_inc_init
190      !!------------------------------------------------------------------------
191
192#if ! defined key_top || ( ! defined key_hadocc && ! defined key_medusa )
193      CALL ctl_stop( ' Attempting to assimilate biogeochemical observations', &
194         &           ' but no compatible biogeochemical model is available' )
195#endif
196
197#if defined key_hadocc
198      IF ( ln_slchldiainc .OR. ln_slchlnoninc .OR. ln_slphydiainc .OR. &
199         & ln_slphynoninc .OR. ln_psi4inc .OR. ln_pphinc .OR. ln_po2inc ) THEN
200         CALL ctl_stop( ' Cannot assimilate PFTs, Si4, pH or O2 into HadOCC' )
201      ENDIF
202#endif
203
204#if defined key_medusa
205      IF ( .NOT. ln_foam_medusa ) THEN
206         CALL ctl_stop( ' MEDUSA data assimilation options not turned on: set ln_foam_medusa' )
207      ENDIF
208#endif
209
210      IF ( ( ln_phytobal ).AND.                                       &
211         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
212         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
213         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
214         & ( .NOT. ln_slphynoninc ) ) THEN
215         CALL ctl_warn( ' Cannot calculate phytoplankton balancing increments', &
216            &           ' if not assimilating ocean colour,',                   &
217            &           ' so ln_phytobal will be set to .false.')
218         ln_phytobal = .FALSE.
219      ENDIF
220
221      IF ( ( ln_balwri ).AND.                                         &
222         & ( .NOT. ln_slchltotinc ).AND.( .NOT. ln_slchldiainc ).AND. &
223         & ( .NOT. ln_slchlnoninc ).AND.( .NOT. ln_schltotinc  ).AND. &
224         & ( .NOT. ln_slphytotinc ).AND.( .NOT. ln_slphydiainc ).AND. &
225         & ( .NOT. ln_slphynoninc ).AND.( .NOT. ln_plchltotinc ).AND. &
226         & ( .NOT. ln_pchltotinc  ).AND.( .NOT. ln_pphinc      ).AND. &
227         & ( .NOT. ln_spco2inc    ).AND.( .NOT. ln_sfco2inc    ) ) THEN
228         CALL ctl_warn( ' Balancing increments not being calculated', &
229            &           ' so ln_balwri will be set to .false.')
230         ln_balwri = .FALSE.
231      ENDIF
232
233      IF ( ln_spco2inc .AND. ln_sfco2inc ) THEN
234         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
235      ENDIF
236
237      IF ( ln_slchltotinc .AND. ln_schltotinc ) THEN
238         CALL ctl_stop( ' Can only assimilate surface log10(chlorophyll) or chlorophyll, not both' )
239      ENDIF
240
241      IF ( ln_plchltotinc .AND. ln_pchltotinc ) THEN
242         CALL ctl_stop( ' Can only assimilate profile log10(chlorophyll) or chlorophyll, not both' )
243      ENDIF
244
245      IF ( ( ln_slchltotinc .OR. ln_schltotinc  ) .AND. &
246         & ( ln_slchldiainc .OR. ln_slchlnoninc ) ) THEN
247         CALL ctl_stop( ' Can only assimilate total or PFT surface chlorophyll, not both' )
248      ENDIF
249
250      IF ( ln_phytobal .AND.                                      &
251         & ( ( ln_slchlnoninc .AND.( .NOT. ln_slchldiainc ) ).OR. &
252         &   ( ln_slchldiainc .AND.( .NOT. ln_slchlnoninc ) ) ) ) THEN
253         CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
254            &           ' unless assimilating all model PFTs,')
255      ENDIF
256
257      IF ( ln_slphytotinc .AND. ( ln_slphydiainc .OR. ln_slphynoninc ) ) THEN
258         CALL ctl_stop( ' Can only assimilate total or PFT surface phytoplankton carbon, not both' )
259      ENDIF
260
261      IF ( ln_phytobal .AND.                                      &
262         & ( ( ln_slphynoninc .AND.( .NOT. ln_slphydiainc ) ).OR. &
263         &   ( ln_slphydiainc .AND.( .NOT. ln_slphynoninc ) ) ) ) THEN
264         CALL ctl_stop( ' Cannot calculate phytoplankton balancing increments', &
265            &           ' unless assimilating all model PFTs,')
266      ENDIF
267
268   END SUBROUTINE asm_bgc_check_options
269
270   !!===========================================================================
271   !!===========================================================================
272   !!===========================================================================
273
274   SUBROUTINE asm_bgc_init_incs( knum )
275      !!------------------------------------------------------------------------
276      !!                    ***  ROUTINE asm_bgc_init_incs  ***
277      !!
278      !! ** Purpose :   initialise bgc increments
279      !!
280      !! ** Method  :   allocate arrays and read increments from file
281      !!
282      !! ** Action  :   allocate arrays and read increments from file
283      !!
284      !! References :   asm_inc_init
285      !!------------------------------------------------------------------------
286      !!
287      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
288      !!
289      !!------------------------------------------------------------------------
290
291      ! Allocate and read increments
292     
293      IF ( ln_slchltotinc ) THEN
294         ALLOCATE( slchltot_bkginc(jpi,jpj) )
295         CALL asm_bgc_read_incs_2d( knum, 'bckinslchltot', slchltot_bkginc )
296      ENDIF
297     
298      IF ( ln_slchldiainc ) THEN
299         ALLOCATE( slchldia_bkginc(jpi,jpj) )
300         CALL asm_bgc_read_incs_2d( knum, 'bckinslchldia', slchldia_bkginc )
301      ENDIF
302     
303      IF ( ln_slchlnoninc ) THEN
304         ALLOCATE( slchlnon_bkginc(jpi,jpj) )
305         CALL asm_bgc_read_incs_2d( knum, 'bckinslchlnon', slchlnon_bkginc )
306      ENDIF
307     
308      IF ( ln_schltotinc ) THEN
309         ALLOCATE( schltot_bkginc(jpi,jpj) )
310         CALL asm_bgc_read_incs_2d( knum, 'bckinschltot', schltot_bkginc )
311      ENDIF
312     
313      IF ( ln_slphytotinc ) THEN
314         ALLOCATE( slphytot_bkginc(jpi,jpj) )
315         CALL asm_bgc_read_incs_2d( knum, 'bckinslphytot', slphytot_bkginc )
316      ENDIF
317     
318      IF ( ln_slphydiainc ) THEN
319         ALLOCATE( slphydia_bkginc(jpi,jpj) )
320         CALL asm_bgc_read_incs_2d( knum, 'bckinslphydia', slphydia_bkginc )
321      ENDIF
322     
323      IF ( ln_slphynoninc ) THEN
324         ALLOCATE( slphynon_bkginc(jpi,jpj) )
325         CALL asm_bgc_read_incs_2d( knum, 'bckinslphynon', slphynon_bkginc )
326      ENDIF
327
328      IF ( ln_sfco2inc ) THEN
329         ALLOCATE( sfco2_bkginc(jpi,jpj) )
330         CALL asm_bgc_read_incs_2d( knum, 'bckinsfco2', sfco2_bkginc )
331      ENDIF
332
333      IF ( ln_spco2inc ) THEN
334         ALLOCATE( spco2_bkginc(jpi,jpj) )
335         CALL asm_bgc_read_incs_2d( knum, 'bckinspco2', sfco2_bkginc )
336      ENDIF
337     
338      IF ( ln_plchltotinc ) THEN
339         ALLOCATE( plchltot_bkginc(jpi,jpj,jpk) )
340         CALL asm_bgc_read_incs_3d( knum, 'bckinplchltot', plchltot_bkginc )
341      ENDIF
342     
343      IF ( ln_pchltotinc ) THEN
344         ALLOCATE( pchltot_bkginc(jpi,jpj,jpk) )
345         CALL asm_bgc_read_incs_3d( knum, 'bckinpchltot', pchltot_bkginc )
346      ENDIF
347     
348      IF ( ln_pno3inc ) THEN
349         ALLOCATE( pno3_bkginc(jpi,jpj,jpk) )
350         CALL asm_bgc_read_incs_3d( knum, 'bckinpno3', pno3_bkginc )
351      ENDIF
352     
353      IF ( ln_psi4inc ) THEN
354         ALLOCATE( psi4_bkginc(jpi,jpj,jpk) )
355         CALL asm_bgc_read_incs_3d( knum, 'bckinpsi4', psi4_bkginc )
356      ENDIF
357     
358      IF ( ln_pdicinc ) THEN
359         ALLOCATE( pdic_bkginc(jpi,jpj,jpk) )
360         CALL asm_bgc_read_incs_3d( knum, 'bckinpdic', pdic_bkginc )
361      ENDIF
362     
363      IF ( ln_palkinc ) THEN
364         ALLOCATE( palk_bkginc(jpi,jpj,jpk) )
365         CALL asm_bgc_read_incs_3d( knum, 'bckinpalk', palk_bkginc )
366      ENDIF
367     
368      IF ( ln_pphinc ) THEN
369         ALLOCATE( pph_bkginc(jpi,jpj,jpk) )
370         CALL asm_bgc_read_incs_3d( knum, 'bckinpph', pph_bkginc )
371      ENDIF
372     
373      IF ( ln_po2inc ) THEN
374         ALLOCATE( po2_bkginc(jpi,jpj,jpk) )
375         CALL asm_bgc_read_incs_3d( knum, 'bckinpo2', po2_bkginc )
376      ENDIF
377
378      ! Allocate balancing increments
379     
380      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
381         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
382         & ln_slphynoninc ) THEN
383#if defined key_top
384         ALLOCATE( phyto2d_balinc(jpi,jpj,jpk,jptra) )
385         phyto2d_balinc(:,:,:,:) = 0.0
386#else
387         CALL ctl_stop( ' key_top must be set for balancing increments' )
388#endif
389      ENDIF
390     
391      IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
392#if defined key_top
393         ALLOCATE( phyto3d_balinc(jpi,jpj,jpk,jptra) )
394         phyto3d_balinc(:,:,:,:) = 0.0
395#else
396         CALL ctl_stop( ' key_top must be set for balancing increments' )
397#endif
398      ENDIF
399
400      IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
401#if defined key_top
402         ALLOCATE( pco2_balinc(jpi,jpj,jpk,jptra) )
403         pco2_balinc(:,:,:,:) = 0.0
404#else
405         CALL ctl_stop( ' key_top must be set for balancing increments' )
406#endif
407      ENDIF
408
409      IF ( ln_pphinc ) THEN
410#if defined key_top
411         ALLOCATE( ph_balinc(jpi,jpj,jpk,jptra) )
412         ph_balinc(:,:,:,:) = 0.0
413#else
414         CALL ctl_stop( ' key_top must be set for balancing increments' )
415#endif
416      ENDIF
417
418   END SUBROUTINE asm_bgc_init_incs
419
420   !!===========================================================================
421   !!===========================================================================
422   !!===========================================================================
423
424   SUBROUTINE asm_bgc_read_incs_2d( knum, cd_bgcname, p_incs )
425      !!------------------------------------------------------------------------
426      !!                    ***  ROUTINE asm_bgc_init_incs  ***
427      !!
428      !! ** Purpose :   read 2d bgc increments
429      !!
430      !! ** Method  :   read increments from file
431      !!
432      !! ** Action  :   read increments from file
433      !!
434      !! References :   asm_inc_init
435      !!------------------------------------------------------------------------
436      !!
437      INTEGER,                      INTENT(in   ) :: knum       ! i/o unit
438      CHARACTER(LEN=13),            INTENT(in   ) :: cd_bgcname ! variable
439      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) :: p_incs     ! increments
440      !!
441      !!------------------------------------------------------------------------
442
443      ! Initialise
444      p_incs(:,:) = 0.0
445     
446      ! read from file
447      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:), 1 )
448     
449      ! Apply the masks
450      p_incs(:,:) = p_incs(:,:) * tmask(:,:,1)
451     
452      ! Set missing increments to 0.0 rather than 1e+20
453      ! to allow for differences in masks
454      WHERE( ABS( p_incs(:,:) ) > 1.0e+10 ) p_incs(:,:) = 0.0
455
456   END SUBROUTINE asm_bgc_read_incs_2d
457
458   !!===========================================================================
459   !!===========================================================================
460   !!===========================================================================
461
462   SUBROUTINE asm_bgc_read_incs_3d( knum, cd_bgcname, p_incs )
463      !!------------------------------------------------------------------------
464      !!                    ***  ROUTINE asm_bgc_init_incs  ***
465      !!
466      !! ** Purpose :   read 3d bgc increments
467      !!
468      !! ** Method  :   read increments from file
469      !!
470      !! ** Action  :   read increments from file
471      !!
472      !! References :   asm_inc_init
473      !!------------------------------------------------------------------------
474      !!
475      INTEGER,                          INTENT(in   ) :: knum       ! i/o unit
476      CHARACTER(LEN=13),                INTENT(in   ) :: cd_bgcname ! variable
477      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) :: p_incs     ! increments
478      !!
479      !!------------------------------------------------------------------------
480
481      ! Initialise
482      p_incs(:,:,:) = 0.0
483     
484      ! read from file
485      CALL iom_get( knum, jpdom_autoglo, TRIM(cd_bgcname), p_incs(:,:,:), 1 )
486     
487      ! Apply the masks
488      p_incs(:,:,:) = p_incs(:,:,:) * tmask(:,:,:)
489     
490      ! Set missing increments to 0.0 rather than 1e+20
491      ! to allow for differences in masks
492      WHERE( ABS( p_incs(:,:,:) ) > 1.0e+10 ) p_incs(:,:,:) = 0.0
493
494   END SUBROUTINE asm_bgc_read_incs_3d
495
496   !!===========================================================================
497   !!===========================================================================
498   !!===========================================================================
499
500   SUBROUTINE asm_bgc_init_bkg
501      !!------------------------------------------------------------------------
502      !!                    ***  ROUTINE asm_bgc_init_bkg  ***
503      !!
504      !! ** Purpose :   initialise bgc background
505      !!
506      !! ** Method  :   allocate arrays and read background from file
507      !!
508      !! ** Action  :   allocate arrays and read background from file
509      !!
510      !! References :   asm_inc_init
511      !!------------------------------------------------------------------------
512      !!
513      INTEGER :: inum   ! i/o unit of background file
514      INTEGER :: jt     ! loop counter
515      !!
516      !!------------------------------------------------------------------------
517
518#if defined key_top && ( defined key_hadocc || defined key_medusa )
519      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
520         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
521         & ln_slphynoninc .OR. ln_plchltotinc .OR. ln_pchltotinc ) THEN
522
523         ALLOCATE( pgrow_avg_bkg(jpi,jpj)        )
524         ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
525         ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
526         ALLOCATE( mld_max_bkg(jpi,jpj)          )
527         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
528         pgrow_avg_bkg(:,:)  = 0.0
529         ploss_avg_bkg(:,:)  = 0.0
530         phyt_avg_bkg(:,:)   = 0.0
531         mld_max_bkg(:,:)    = 0.0
532         tracer_bkg(:,:,:,:) = 0.0
533
534#if defined key_hadocc
535         ALLOCATE( chl_bkg(jpi,jpj,jpk)    )
536         ALLOCATE( cchl_p_bkg(jpi,jpj,jpk) )
537         chl_bkg(:,:,:)    = 0.0
538         cchl_p_bkg(:,:,:) = 0.0
539#endif
540         
541         !--------------------------------------------------------------------
542         ! Read background variables for phytoplankton assimilation
543         ! Some only required if performing balancing
544         !--------------------------------------------------------------------
545
546         CALL iom_open( c_asmbkg, inum )
547
548#if defined key_hadocc
549         CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl',  chl_bkg    )
550         CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
551         chl_bkg(:,:,:)    = chl_bkg(:,:,:)    * tmask(:,:,:)
552         cchl_p_bkg(:,:,:) = cchl_p_bkg(:,:,:) * tmask(:,:,:)
553#elif defined key_medusa
554         CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
555         CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
556         CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
557         CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
558         CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
559#endif
560         
561         IF ( ln_phytobal ) THEN
562
563            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
564            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
565            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  )
566            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
567            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
568            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
569            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1)
570            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1)
571
572#if defined key_hadocc
573            CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
574            CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
575            CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
576            CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
577            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
578            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
579#elif defined key_medusa
580            CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
581            CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
582            CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
583            CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
584            CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
585            CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
586            CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
587            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
588            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
589            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
590#endif
591         ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
592#if defined key_hadocc
593            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
594            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
595#elif defined key_medusa
596            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
597            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
598#endif
599            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
600            mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
601         ENDIF
602
603         CALL iom_close( inum )
604         
605         DO jt = 1, jptra
606            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
607         END DO
608     
609      ELSE IF ( ln_spco2inc .OR. ln_sfco2inc .OR. ln_pphinc ) THEN
610
611         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
612         ALLOCATE( mld_max_bkg(jpi,jpj)          )
613         tracer_bkg(:,:,:,:) = 0.0
614         mld_max_bkg(:,:)    = 0.0
615
616         CALL iom_open( c_asmbkg, inum )
617         
618#if defined key_hadocc
619         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
620         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
621#elif defined key_medusa
622         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
623         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
624#endif
625         CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
626
627         CALL iom_close( inum )
628         
629         DO jt = 1, jptra
630            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
631         END DO
632         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
633     
634      ENDIF
635#else
636      CALL ctl_stop( ' asm_bgc_init_bkg: key_top and a compatible biogeochemical model required' )
637#endif
638
639   END SUBROUTINE asm_bgc_init_bkg
640
641   !!===========================================================================
642   !!===========================================================================
643   !!===========================================================================
644
645   SUBROUTINE asm_bgc_bal_wri( kt )
646      !!------------------------------------------------------------------------
647      !!
648      !!                  ***  ROUTINE asm_bgc_bal_wri ***
649      !!
650      !! ** Purpose : Write to file the balancing increments
651      !!              calculated for biogeochemistry
652      !!
653      !! ** Method  : Write to file the balancing increments
654      !!              calculated for biogeochemistry
655      !!
656      !! ** Action  :
657      !!                   
658      !! References :
659      !!
660      !! History    :
661      !!        ! 2014-08 (D. Ford)  Initial version, based on asm_bkg_wri
662      !!------------------------------------------------------------------------
663      !! * Arguments
664      INTEGER, INTENT( IN ) :: kt        ! Current time-step
665      !! * Local declarations
666      CHARACTER(LEN=50) :: cl_asmbal     ! Filename (with extension)
667      LOGICAL           :: llok          ! Check if file exists
668      INTEGER           :: inum          ! File unit number
669      REAL(wp)          :: zdate         ! Date
670      !!------------------------------------------------------------------------
671     
672      ! Set things up
673      zdate = REAL( ndastp )
674      WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
675
676      ! Check if file exists
677      INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
678      IF( .NOT. llok ) THEN
679         IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
680            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
681
682         ! Define the output file       
683         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
684
685         ! Write the information
686         CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   )
687
688         IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
689            & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
690            & ln_slphynoninc ) THEN
691#if defined key_medusa
692            CALL iom_rstput( kt, kt, inum, 'phy2d_chn', phyto2d_balinc(:,:,:,jpchn) )
693            CALL iom_rstput( kt, kt, inum, 'phy2d_chd', phyto2d_balinc(:,:,:,jpchd) )
694            CALL iom_rstput( kt, kt, inum, 'phy2d_phn', phyto2d_balinc(:,:,:,jpphn) )
695            CALL iom_rstput( kt, kt, inum, 'phy2d_phd', phyto2d_balinc(:,:,:,jpphd) )
696            CALL iom_rstput( kt, kt, inum, 'phy2d_pds', phyto2d_balinc(:,:,:,jppds) )
697            IF ( ln_phytobal ) THEN
698               CALL iom_rstput( kt, kt, inum, 'phy2d_zmi', phyto2d_balinc(:,:,:,jpzmi) )
699               CALL iom_rstput( kt, kt, inum, 'phy2d_zme', phyto2d_balinc(:,:,:,jpzme) )
700               CALL iom_rstput( kt, kt, inum, 'phy2d_din', phyto2d_balinc(:,:,:,jpdin) )
701               CALL iom_rstput( kt, kt, inum, 'phy2d_sil', phyto2d_balinc(:,:,:,jpsil) )
702               CALL iom_rstput( kt, kt, inum, 'phy2d_fer', phyto2d_balinc(:,:,:,jpfer) )
703               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jpdet) )
704               CALL iom_rstput( kt, kt, inum, 'phy2d_dtc', phyto2d_balinc(:,:,:,jpdtc) )
705               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jpdic) )
706               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jpalk) )
707               CALL iom_rstput( kt, kt, inum, 'phy2d_oxy', phyto2d_balinc(:,:,:,jpoxy) )
708            ENDIF
709#elif defined key_hadocc
710            CALL iom_rstput( kt, kt, inum, 'phy2d_phy', phyto2d_balinc(:,:,:,jp_had_phy) )
711            IF ( ln_phytobal ) THEN
712               CALL iom_rstput( kt, kt, inum, 'phy2d_nut', phyto2d_balinc(:,:,:,jp_had_nut) )
713               CALL iom_rstput( kt, kt, inum, 'phy2d_zoo', phyto2d_balinc(:,:,:,jp_had_zoo) )
714               CALL iom_rstput( kt, kt, inum, 'phy2d_det', phyto2d_balinc(:,:,:,jp_had_pdn) )
715               CALL iom_rstput( kt, kt, inum, 'phy2d_dic', phyto2d_balinc(:,:,:,jp_had_dic) )
716               CALL iom_rstput( kt, kt, inum, 'phy2d_alk', phyto2d_balinc(:,:,:,jp_had_alk) )
717            ENDIF
718#endif
719         ENDIF
720
721         IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
722#if defined key_medusa
723            CALL iom_rstput( kt, kt, inum, 'phy3d_chn', phyto3d_balinc(:,:,:,jpchn) )
724            CALL iom_rstput( kt, kt, inum, 'phy3d_chd', phyto3d_balinc(:,:,:,jpchd) )
725            CALL iom_rstput( kt, kt, inum, 'phy3d_phn', phyto3d_balinc(:,:,:,jpphn) )
726            CALL iom_rstput( kt, kt, inum, 'phy3d_phd', phyto3d_balinc(:,:,:,jpphd) )
727            CALL iom_rstput( kt, kt, inum, 'phy3d_pds', phyto3d_balinc(:,:,:,jppds) )
728#elif defined key_hadocc
729            CALL iom_rstput( kt, kt, inum, 'phy3d_phy', phyto3d_balinc(:,:,:,jp_had_phy) )
730#endif
731         ENDIF
732
733         IF ( ln_spco2inc ) THEN
734#if defined key_medusa
735            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jpdic) )
736            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jpalk) )
737#elif defined key_hadocc
738            CALL iom_rstput( kt, kt, inum, 'pco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
739            CALL iom_rstput( kt, kt, inum, 'pco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
740#endif
741         ELSE IF ( ln_sfco2inc ) THEN
742#if defined key_medusa
743            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jpdic) )
744            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jpalk) )
745#elif defined key_hadocc
746            CALL iom_rstput( kt, kt, inum, 'fco2_dic', pco2_balinc(:,:,:,jp_had_dic) )
747            CALL iom_rstput( kt, kt, inum, 'fco2_alk', pco2_balinc(:,:,:,jp_had_alk) )
748#endif
749         ENDIF
750
751         IF ( ln_pphinc ) THEN
752#if defined key_medusa
753            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jpdic) )
754            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jpalk) )
755#elif defined key_hadocc
756            CALL iom_rstput( kt, kt, inum, 'ph_dic', ph_balinc(:,:,:,jp_had_dic) )
757            CALL iom_rstput( kt, kt, inum, 'ph_alk', ph_balinc(:,:,:,jp_had_alk) )
758#endif
759         ENDIF
760
761         CALL iom_close( inum )
762      ELSE
763         CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
764            &           ' Therefore not writing out balancing increments at this timestep', &
765            &           ' Something has probably gone wrong somewhere' )
766         IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
767            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
768      ENDIF
769                                 
770   END SUBROUTINE asm_bgc_bal_wri
771
772   !!===========================================================================
773   !!===========================================================================
774   !!===========================================================================
775
776   SUBROUTINE asm_bgc_bkg_wri( kt, knum )
777      !!------------------------------------------------------------------------
778      !!                    ***  ROUTINE asm_bgc_bkg_wri  ***
779      !!
780      !! ** Purpose :   write out bgc background
781      !!
782      !! ** Method  :   write out bgc background
783      !!
784      !! ** Action  :   write out bgc background
785      !!
786      !! References :   asm_bkg_wri
787      !!------------------------------------------------------------------------
788      !!
789      INTEGER, INTENT(in   ) :: kt     ! Current time-step
790      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
791      !!
792      !!------------------------------------------------------------------------
793
794#if defined key_hadocc
795      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg             )
796      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg             )
797      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg              )
798      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max               )
799      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) )
800      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) )
801      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) )
802      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) )
803      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) )
804      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) )
805      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,:)     )
806      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         )
807#elif defined key_medusa
808      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg        )
809      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        )
810      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         )
811      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          )
812      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) )
813      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd'  , trn(:,:,:,jpchd) )
814      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn'  , trn(:,:,:,jpphn) )
815      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd'  , trn(:,:,:,jpphd) )
816      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds'  , trn(:,:,:,jppds) )
817      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi'  , trn(:,:,:,jpzmi) )
818      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme'  , trn(:,:,:,jpzme) )
819      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din'  , trn(:,:,:,jpdin) )
820      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil'  , trn(:,:,:,jpsil) )
821      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer'  , trn(:,:,:,jpfer) )
822      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det'  , trn(:,:,:,jpdet) )
823      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc'  , trn(:,:,:,jpdtc) )
824      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic'  , trn(:,:,:,jpdic) )
825      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk'  , trn(:,:,:,jpalk) )
826      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
827#endif
828
829   END SUBROUTINE asm_bgc_bkg_wri
830
831   !!===========================================================================
832   !!===========================================================================
833   !!===========================================================================
834
835   SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog )
836      !!------------------------------------------------------------------------
837      !!                    ***  ROUTINE asm_bgc_init_incs  ***
838      !!
839      !! ** Purpose :   convert log increments to non-log
840      !!
841      !! ** Method  :   need to account for model background,
842      !!                cannot simply do 10^log_inc. Need to:
843      !!                1) Add log_inc to log10(background) to get log10(analysis)
844      !!                2) Take 10^log10(analysis) to get analysis
845      !!                3) Subtract background from analysis to get non-log incs
846      !!
847      !! ** Action  :   return non-log increments
848      !!
849      !! References :   
850      !!------------------------------------------------------------------------
851      !!
852      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pbkg        ! Background
853      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pinc_log    ! Log incs
854      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs
855      !
856      INTEGER                                     :: ji, jj      ! Loop counters
857      !!
858      !!------------------------------------------------------------------------
859
860      DO jj = 1, jpj
861         DO ji = 1, jpi
862            IF ( pbkg(ji,jj) > 0.0 ) THEN
863               pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + &
864                  &                       pinc_log(ji,jj) )      &
865                  &                 - pbkg(ji,jj)
866            ELSE
867               pinc_nonlog(ji,jj) = 0.0
868            ENDIF
869         END DO
870      END DO
871
872   END SUBROUTINE asm_bgc_unlog_2d
873
874   !!===========================================================================
875   !!===========================================================================
876   !!===========================================================================
877
878   SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
879      !!------------------------------------------------------------------------
880      !!                    ***  ROUTINE phyto2d_asm_inc  ***
881      !!         
882      !! ** Purpose : Apply the chlorophyll assimilation increments.
883      !!
884      !! ** Method  : Calculate increments to state variables using nitrogen
885      !!              balancing.
886      !!              Direct initialization or Incremental Analysis Updating.
887      !!
888      !! ** Action  :
889      !!------------------------------------------------------------------------
890      INTEGER,  INTENT(IN) :: kt        ! Current time step
891      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
892      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
893      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
894      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
895      !
896      INTEGER                      :: jk          ! Loop counter
897      INTEGER                      :: it          ! Index
898      REAL(wp)                     :: zincwgt     ! IAU weight for current time step
899      REAL(wp)                     :: zincper     ! IAU interval in seconds
900      REAL(wp), DIMENSION(jpi,jpj) :: zmld        ! Mixed layer depth
901      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs
902      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
903      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs
904      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
905#if defined key_medusa
906      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
907      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
908      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs
909      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
910      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs
911      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
912      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
913      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
914#endif
915      !!------------------------------------------------------------------------
916     
917      IF ( kt <= nit000 ) THEN
918
919         ! Un-log any log increments for passing to balancing routines
920         ! Remember that two sets of non-log increments should not be
921         ! expected to be in the same ratio as their log equivalents
922         
923         ! Total chlorophyll
924         IF ( ln_slchltotinc ) THEN
925#if defined key_medusa
926            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
927#elif defined key_hadocc
928            zbkg_chltot(:,:) = chl_bkg(:,:,1)
929#endif
930            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
931         ELSE IF ( ln_schltotinc ) THEN
932            zinc_chltot(:,:) = schltot_bkginc(:,:)
933         ELSE
934            zinc_chltot(:,:) = 0.0
935         ENDIF
936
937#if defined key_medusa
938         ! Diatom chlorophyll
939         IF ( ln_slchldiainc ) THEN
940            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
941            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
942         ELSE
943            zinc_chldia(:,:) = 0.0
944         ENDIF
945#endif
946
947#if defined key_medusa
948         ! Non-diatom chlorophyll
949         IF ( ln_slchlnoninc ) THEN
950            zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
951            CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon )
952         ELSE
953            zinc_chlnon(:,:) = 0.0
954         ENDIF
955#endif
956
957         ! Total phytoplankton carbon
958         IF ( ln_slphytotinc ) THEN
959#if defined key_medusa
960            zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
961#elif defined key_hadocc
962            zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
963#endif
964            CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot )
965         ELSE
966            zinc_phytot(:,:) = 0.0
967         ENDIF
968
969#if defined key_medusa
970         ! Diatom phytoplankton carbon
971         IF ( ln_slphydiainc ) THEN
972            zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
973            CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia )
974         ELSE
975            zinc_phydia(:,:) = 0.0
976         ENDIF
977#endif
978
979#if defined key_medusa
980         ! Non-diatom phytoplankton carbon
981         IF ( ln_slphynoninc ) THEN
982            zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
983            CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon )
984         ELSE
985            zinc_phynon(:,:) = 0.0
986         ENDIF
987#endif
988
989         ! Select mixed layer
990         IF ( ll_asmdin ) THEN
991#if defined key_top && ( defined key_hadocc || defined key_medusa )
992            CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
993               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
994            zmld(:,:) = mld_max_bkg(:,:)
995#else
996            CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
997#endif
998         ELSE
999            SELECT CASE( mld_choice_bgc )
1000            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
1001               zmld(:,:) = hmld(:,:)
1002            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1003               zmld(:,:) = hmlp(:,:)
1004            CASE ( 3 )                   ! Kara MLD [Interpolated]
1005#if defined key_karaml
1006               IF ( ln_kara ) THEN
1007                  zmld(:,:) = hmld_kara(:,:)
1008               ELSE
1009                  CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1010                     &           ' but ln_kara=.false.' )
1011               ENDIF
1012#else
1013               CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1014                  &           ' but is not defined' )
1015#endif
1016            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
1017               !zmld(:,:) = hmld_tref(:,:)
1018               CALL ctl_stop( ' hmld_tref mixed layer requested for phyto2d assimilation,', &
1019                  &           ' but is not available in this version' )
1020            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1021               zmld(:,:) = hmlpt(:,:)
1022            END SELECT
1023         ENDIF
1024
1025         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
1026
1027#if defined key_medusa
1028         CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
1029            &                         zinc_chltot,                         &
1030            &                         ln_slchldiainc,                      &
1031            &                         zinc_chldia,                         &
1032            &                         ln_slchlnoninc,                      &
1033            &                         zinc_chlnon,                         &
1034            &                         ln_slphytotinc,                      &
1035            &                         zinc_phytot,                         &
1036            &                         ln_slphydiainc,                      &
1037            &                         zinc_phydia,                         &
1038            &                         ln_slphynoninc,                      &
1039            &                         zinc_phynon,                         &
1040            &                         zincper,                             &
1041            &                         rn_maxchlinc, ln_phytobal, zmld,     &
1042            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
1043            &                         phyt_avg_bkg, mld_max_bkg,           &
1044            &                         tracer_bkg, phyto2d_balinc )
1045#elif defined key_hadocc
1046         CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
1047            &                         zinc_chltot,                         &
1048            &                         ln_slphytotinc,                      &
1049            &                         zinc_phytot,                         &
1050            &                         zincper,                             &
1051            &                         rn_maxchlinc, ln_phytobal, zmld,     &
1052            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
1053            &                         phyt_avg_bkg, mld_max_bkg,           &
1054            &                         cchl_p_bkg(:,:,1),                   &
1055            &                         tracer_bkg, phyto2d_balinc )
1056#else
1057         CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
1058            &           'but not defined a biogeochemical model' )
1059#endif
1060
1061      ENDIF
1062
1063      IF ( ll_asmiau ) THEN
1064
1065         !--------------------------------------------------------------------
1066         ! Incremental Analysis Updating
1067         !--------------------------------------------------------------------
1068
1069         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1070
1071            it = kt - nit000 + 1
1072            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1073            ! note this is not a tendency so should not be divided by rdt
1074
1075            IF(lwp) THEN
1076               WRITE(numout,*) 
1077               WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
1078                  &  kt,' with IAU weight = ', pwgtiau(it)
1079               WRITE(numout,*) '~~~~~~~~~~~~'
1080            ENDIF
1081
1082            ! Update the biogeochemical variables
1083            ! Add directly to trn and trb, rather than to tra, because tra gets
1084            ! reset to zero at the start of trc_stp, called after this routine
1085#if defined key_medusa
1086            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1087               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1088               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1089                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1090               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1091                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1092            END WHERE
1093#elif defined key_hadocc
1094            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1095               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1096               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1097                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1098               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1099                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1100            END WHERE
1101#endif
1102
1103            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1104            ! which is called at end of model run
1105         ENDIF
1106
1107      ELSEIF ( ll_asmdin ) THEN 
1108
1109         !--------------------------------------------------------------------
1110         ! Direct Initialization
1111         !--------------------------------------------------------------------
1112         
1113         IF ( kt == nitdin_r ) THEN
1114
1115            neuler = 0                    ! Force Euler forward step
1116
1117            ! Initialize the now fields with the background + increment
1118            ! Background currently is what the model is initialised with
1119            CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
1120               &           ' Background state is taken from model rather than background file' )
1121#if defined key_medusa
1122            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1123               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1124               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1125                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
1126               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1127            END WHERE
1128#elif defined key_hadocc
1129            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1130               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1131               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1132                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1)
1133               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1134            END WHERE
1135#endif
1136 
1137            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1138            ! which is called at end of model run
1139         ENDIF
1140         !
1141      ENDIF
1142      !
1143   END SUBROUTINE phyto2d_asm_inc
1144
1145   !!===========================================================================
1146   !!===========================================================================
1147   !!===========================================================================
1148
1149   SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1150      !!------------------------------------------------------------------------
1151      !!                    ***  ROUTINE phyto3d_asm_inc  ***
1152      !!         
1153      !! ** Purpose : Apply the profile chlorophyll assimilation increments.
1154      !!
1155      !! ** Method  : Calculate increments to state variables.
1156      !!              Direct initialization or Incremental Analysis Updating.
1157      !!
1158      !! ** Action  :
1159      !!------------------------------------------------------------------------
1160      INTEGER,  INTENT(IN) :: kt        ! Current time step
1161      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1162      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1163      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1164      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1165      !
1166      INTEGER                          :: ji, jj, jk   ! Loop counters
1167      INTEGER                          :: it           ! Index
1168      REAL(wp)                         :: zincwgt      ! IAU weight for timestep
1169      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn
1170      REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd
1171      REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio
1172      REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio
1173      REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio
1174      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments
1175      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background
1176      !!------------------------------------------------------------------------
1177
1178      IF ( kt <= nit000 ) THEN
1179
1180         IF ( ln_plchltotinc ) THEN
1181            ! Convert log10(chlorophyll) increment back to a chlorophyll increment
1182            ! In order to transform logchl incs to chl incs, need to account for model
1183            ! background, cannot simply do 10^logchl_bkginc. Need to:
1184            ! 1) Add logchl inc to log10(background) to get log10(analysis)
1185            ! 2) Take 10^log10(analysis) to get analysis
1186            ! 3) Subtract background from analysis to get chl incs
1187            ! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
1188#if defined key_medusa
1189            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
1190#elif defined key_hadocc
1191            bkg_chl(:,:,:) = chl_bkg(:,:,:)
1192#endif
1193            DO jk = 1, jpk
1194               DO jj = 1, jpj
1195                  DO ji = 1, jpi
1196                     IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
1197                        chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
1198                        IF ( rn_maxchlinc > 0.0 ) THEN
1199                           chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
1200                        ENDIF
1201                     ELSE
1202                        chl_inc(ji,jj,jk) = 0.0
1203                     ENDIF
1204                  END DO
1205               END DO
1206            END DO
1207         ELSE IF ( ln_pchltotinc ) THEN
1208            DO jk = 1, jpk
1209               DO jj = 1, jpj
1210                  DO ji = 1, jpi
1211                     IF ( rn_maxchlinc > 0.0 ) THEN
1212                        chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
1213                     ELSE
1214                        chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
1215                     ENDIF
1216                  END DO
1217               END DO
1218            END DO
1219         ENDIF
1220
1221#if defined key_medusa
1222         ! Loop over each grid point partioning the increments based on existing ratios
1223         DO jk = 1, jpk
1224            DO jj = 1, jpj
1225               DO ji = 1, jpi
1226                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
1227                     zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
1228                     zfrac_chd = 1.0 - zfrac_chn
1229                     phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
1230                     phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
1231                     zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
1232                     zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
1233                     zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
1234                     phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
1235                     phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
1236                     phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
1237                  ENDIF
1238               END DO
1239            END DO
1240         END DO
1241#elif defined key_hadocc
1242         phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
1243#else
1244         CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
1245            &           'but not defined a biogeochemical model' )
1246#endif
1247
1248      ENDIF
1249
1250      IF ( ll_asmiau ) THEN
1251
1252         !--------------------------------------------------------------------
1253         ! Incremental Analysis Updating
1254         !--------------------------------------------------------------------
1255
1256         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1257
1258            it = kt - nit000 + 1
1259            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1260            ! note this is not a tendency so should not be divided by rdt
1261
1262            IF(lwp) THEN
1263               WRITE(numout,*) 
1264               WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
1265                  &  kt,' with IAU weight = ', pwgtiau(it)
1266               WRITE(numout,*) '~~~~~~~~~~~~'
1267            ENDIF
1268
1269            ! Update the biogeochemical variables
1270            ! Add directly to trn and trb, rather than to tra, because tra gets
1271            ! reset to zero at the start of trc_stp, called after this routine
1272#if defined key_medusa
1273            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1274               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1275               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1276                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1277               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1278                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1279            END WHERE
1280#elif defined key_hadocc
1281            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1282               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1283               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1284                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1285               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1286                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1287            END WHERE
1288#endif
1289
1290            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1291            ! which is called at end of model run
1292         ENDIF
1293
1294      ELSEIF ( ll_asmdin ) THEN 
1295
1296         !--------------------------------------------------------------------
1297         ! Direct Initialization
1298         !--------------------------------------------------------------------
1299         
1300         IF ( kt == nitdin_r ) THEN
1301
1302            neuler = 0                    ! Force Euler forward step
1303
1304            ! Initialize the now fields with the background + increment
1305            ! Background currently is what the model is initialised with
1306            CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
1307               &           ' Background state is taken from model rather than background file' )
1308#if defined key_medusa
1309            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1310               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1311               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1312                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
1313               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1314            END WHERE
1315#elif defined key_hadocc
1316            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1317               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1318               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1319                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1)
1320               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1321            END WHERE
1322#endif
1323 
1324            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1325            ! which is called at end of model run
1326         ENDIF
1327         !
1328      ENDIF
1329      !
1330   END SUBROUTINE phyto3d_asm_inc
1331
1332   !!===========================================================================
1333   !!===========================================================================
1334   !!===========================================================================
1335
1336   SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1337      &                     ll_trainc, pt_bkginc, ps_bkginc )
1338      !!------------------------------------------------------------------------
1339      !!                    ***  ROUTINE pco2_asm_inc  ***
1340      !!         
1341      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
1342      !!
1343      !! ** Method  : Calculate increments to state variables using carbon
1344      !!              balancing.
1345      !!              Direct initialization or Incremental Analysis Updating.
1346      !!
1347      !! ** Action  :
1348      !!------------------------------------------------------------------------
1349      INTEGER, INTENT(IN)                   :: kt           ! Current time step
1350      LOGICAL, INTENT(IN)                   :: ll_asmdin    ! Flag for direct initialisation
1351      LOGICAL, INTENT(IN)                   :: ll_asmiau    ! Flag for incremental analysis update
1352      INTEGER, INTENT(IN)                   :: kcycper      ! Dimension of pwgtiau
1353      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1354      LOGICAL, INTENT(IN)                   :: ll_trainc    ! Flag for T/S increments
1355      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1356      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1357      !
1358      INTEGER                               :: ji, jj, jk   ! Loop counters
1359      INTEGER                               :: it           ! Index
1360      INTEGER                               :: jkmax        ! Index of mixed layer depth
1361      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
1362      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
1363      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
1364      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
1365      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
1366      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
1367      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
1368      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
1369      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
1370
1371      ! Coefficients for fCO2 to pCO2 conversion
1372      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
1373      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
1374      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
1375      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
1376      REAL(wp)                              :: zcoef_fco2_5 = 2.0
1377      REAL(wp)                              :: zcoef_fco2_6 = 57.7
1378      REAL(wp)                              :: zcoef_fco2_7 = 0.118
1379      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
1380      !!------------------------------------------------------------------------
1381
1382      IF ( kt <= nit000 ) THEN
1383
1384         !----------------------------------------------------------------------
1385         ! DIC and alkalinity balancing
1386         !----------------------------------------------------------------------
1387
1388         IF ( ln_sfco2inc ) THEN
1389#if defined key_medusa && defined key_roam
1390            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
1391            patm(1) = 1.0
1392            phyd(1) = 0.0
1393            DO jj = 1, jpj
1394               DO ji = 1, jpi
1395                  CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
1396               END DO
1397            END DO
1398#else
1399            ! If assimilating fCO2, then convert to pCO2 using temperature
1400            ! See flux_gas.F90 within HadOCC for details of calculation
1401            pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) /                                                            &
1402               &                    EXP((zcoef_fco2_1                                                            + &
1403               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
1404               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
1405               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
1406               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
1407               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
1408#endif
1409         ELSE
1410            pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
1411         ENDIF
1412
1413         ! Account for physics assimilation if required
1414         IF ( ll_trainc ) THEN
1415            tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
1416            sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
1417         ELSE
1418            tem_bkg_temp(:,:) = tsn(:,:,1,1)
1419            sal_bkg_temp(:,:) = tsn(:,:,1,2)
1420         ENDIF
1421
1422#if defined key_medusa
1423         ! Account for phytoplankton balancing if required
1424         IF ( ln_phytobal ) THEN
1425            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
1426            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
1427         ELSE
1428            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1429            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1430         ENDIF
1431
1432         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1433            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1434            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1435
1436#elif defined key_hadocc
1437         ! Account for phytoplankton balancing if required
1438         IF ( ln_phytobal ) THEN
1439            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
1440            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
1441         ELSE
1442            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1443            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1444         ENDIF
1445
1446         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1447            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1448            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1449
1450#else
1451         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1452            &           'but not defined a biogeochemical model' )
1453#endif
1454
1455         ! Select mixed layer
1456         IF ( ll_asmdin ) THEN
1457#if defined key_hadocc || defined key_medusa
1458            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1459               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1460            zmld(:,:) = mld_max_bkg(:,:)
1461#else
1462            CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
1463#endif
1464         ELSE
1465            SELECT CASE( mld_choice_bgc )
1466            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
1467               zmld(:,:) = hmld(:,:)
1468            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1469               zmld(:,:) = hmlp(:,:)
1470            CASE ( 3 )                   ! Kara MLD [Interpolated]
1471#if defined key_karaml
1472               IF ( ln_kara ) THEN
1473                  zmld(:,:) = hmld_kara(:,:)
1474               ELSE
1475                  CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1476                     &           ' but ln_kara=.false.' )
1477               ENDIF
1478#else
1479               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1480                  &           ' but is not defined' )
1481#endif
1482            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
1483               !zmld(:,:) = hmld_tref(:,:)
1484               CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
1485                  &           ' but is not available in this version' )
1486            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1487               zmld(:,:) = hmlpt(:,:)
1488            END SELECT
1489         ENDIF
1490
1491         ! Propagate through mixed layer
1492         DO jj = 1, jpj
1493            DO ji = 1, jpi
1494               !
1495               jkmax = jpk-1
1496               DO jk = jpk-1, 1, -1
1497                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
1498                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
1499                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
1500                     jkmax = jk
1501                  ENDIF
1502               END DO
1503               !
1504#if defined key_top
1505               DO jk = 2, jkmax
1506                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
1507               END DO
1508#endif
1509               !
1510            END DO
1511         END DO
1512
1513      ENDIF
1514
1515      IF ( ll_asmiau ) THEN
1516
1517         !--------------------------------------------------------------------
1518         ! Incremental Analysis Updating
1519         !--------------------------------------------------------------------
1520
1521         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1522
1523            it = kt - nit000 + 1
1524            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1525            ! note this is not a tendency so should not be divided by rdt
1526
1527            IF(lwp) THEN
1528               WRITE(numout,*) 
1529               IF ( ln_spco2inc ) THEN
1530                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
1531                     &  kt,' with IAU weight = ', pwgtiau(it)
1532               ELSE IF ( ln_sfco2inc ) THEN
1533                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
1534                     &  kt,' with IAU weight = ', pwgtiau(it)
1535               ENDIF
1536               WRITE(numout,*) '~~~~~~~~~~~~'
1537            ENDIF
1538
1539            ! Update the biogeochemical variables
1540            ! Add directly to trn and trb, rather than to tra, because tra gets
1541            ! reset to zero at the start of trc_stp, called after this routine
1542#if defined key_medusa
1543            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1544               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1545               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1546                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1547               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1548                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1549            END WHERE
1550#elif defined key_hadocc
1551            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1552               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1553               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1554                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1555               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1556                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1557            END WHERE
1558#endif
1559
1560            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1561            ! which is called at end of model run
1562
1563         ENDIF
1564
1565      ELSEIF ( ll_asmdin ) THEN 
1566
1567         !--------------------------------------------------------------------
1568         ! Direct Initialization
1569         !--------------------------------------------------------------------
1570         
1571         IF ( kt == nitdin_r ) THEN
1572
1573            neuler = 0                    ! Force Euler forward step
1574
1575            ! Initialize the now fields with the background + increment
1576            ! Background currently is what the model is initialised with
1577            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1578               &           ' Background state is taken from model rather than background file' )
1579#if defined key_medusa
1580            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1581               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1582               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1583                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1584               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1585            END WHERE
1586#elif defined key_hadocc
1587            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1588               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1589               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1590                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
1591               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1592            END WHERE
1593#endif
1594 
1595            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1596            ! which is called at end of model run
1597         ENDIF
1598         !
1599      ENDIF
1600      !
1601   END SUBROUTINE pco2_asm_inc
1602
1603   !!===========================================================================
1604   !!===========================================================================
1605   !!===========================================================================
1606
1607   SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1608      &                   ll_trainc, pt_bkginc, ps_bkginc )
1609      !!------------------------------------------------------------------------
1610      !!                    ***  ROUTINE ph_asm_inc  ***
1611      !!         
1612      !! ** Purpose : Apply the pH assimilation increments.
1613      !!
1614      !! ** Method  : Calculate increments to DIC and alkalinity from pH
1615      !!              Use a similar approach to the pCO2 scheme
1616      !!              Direct initialization or Incremental Analysis Updating.
1617      !!
1618      !! ** Action  :
1619      !!------------------------------------------------------------------------
1620      INTEGER,                          INTENT(IN) :: kt        ! Current time step
1621      LOGICAL,                          INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1622      LOGICAL,                          INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1623      INTEGER,                          INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1624      REAL(wp), DIMENSION(kcycper),     INTENT(IN) :: pwgtiau   ! IAU weights
1625      LOGICAL,                          INTENT(IN) :: ll_trainc ! Flag for T/S increments
1626      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1627      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1628     
1629      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
1630      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
1631      REAL(wp)                         :: PH_OUT         ! pH from pH calculation
1632      REAL(wp)                         :: ph_bkg         ! pH from pH calculation
1633      REAL(wp)                         :: ph_dic         ! pH from pH calculation
1634      REAL(wp)                         :: ph_alk         ! pH from pH calculation
1635      REAL(wp)                         :: dph_ddic       ! Change in pH wrt DIC
1636      REAL(wp)                         :: dph_dalk       ! Change in pH wrt alk
1637      REAL(wp)                         :: weight         ! Increment weighting
1638      REAL(wp)                         :: zincwgt        ! IAU weight for current time step
1639      REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp   ! DIC background
1640      REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp   ! Alkalinity background
1641      REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp   ! DIN background
1642      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp   ! Silicate background
1643      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp   ! Temperature background
1644      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp   ! Salinity background
1645      REAL(wp), DIMENSION(19)          :: dummy_out      ! Dummy variables
1646      INTEGER                          :: ji, jj, jk, jx ! Loop counters
1647      INTEGER                          :: it             ! Index
1648      !!------------------------------------------------------------------------
1649
1650#if ! defined key_medusa
1651      CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
1652#else
1653
1654      IF ( kt <= nit000 ) THEN
1655
1656         !----------------------------------------------------------------------
1657         ! DIC and alkalinity balancing
1658         !----------------------------------------------------------------------
1659
1660         ! Account for physics assimilation if required
1661         IF ( ll_trainc ) THEN
1662            tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
1663            sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
1664         ELSE
1665            tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
1666            sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
1667         ENDIF
1668
1669         ! Account for phytoplankton balancing if required
1670         IF ( ln_phytobal ) THEN
1671            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
1672            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
1673            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
1674            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
1675         ELSE
1676            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
1677            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
1678            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
1679            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
1680         ENDIF
1681         
1682         ! Account for pCO2 balancing if required
1683         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1684            dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
1685            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
1686         ENDIF
1687         
1688         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
1689         ! This requires three calls to the MOCSY carbonate package
1690         ! One to find pH at background DIC and alkalinity
1691         ! One to find pH when DIC is increased by 10
1692         ! One to find pH when alkalinity is increased by 10
1693         ! Then calculate the gradients
1694
1695         DO jk = 1, jpk
1696            DO jj = 1, jpj
1697               DO ji = 1, jpi
1698
1699                  IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
1700
1701                     DO jx = 1, 3
1702                        IF ( jx == 1 ) THEN
1703                           DIC_IN = dic_bkg_temp(ji,jj,jk)
1704                           ALK_IN = alk_bkg_temp(ji,jj,jk)
1705                        ELSE IF ( jx == 2 ) THEN
1706                           DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
1707                           ALK_IN = alk_bkg_temp(ji,jj,jk)
1708                        ELSE IF ( jx == 3 ) THEN
1709                           DIC_IN = dic_bkg_temp(ji,jj,jk)
1710                           ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
1711                        ENDIF
1712
1713                        CALL mocsy_interface(tsn(ji,jj,jk,jp_tem),         & ! IN:  temperature
1714                           &                 tsn(ji,jj,jk,jp_sal),         & !      salinity
1715                           &                 ALK_IN,                       & !      alkalinity
1716                           &                 DIC_IN,                       & !      DIC
1717                           &                 sil_bkg_temp(ji,jj,jk),       & !      silicate
1718                           &                 din_bkg_temp(ji,jj,jk)/16.0,  & !      phosphate (Redfield from DIN)
1719                           &                 1.0,                          & !      atmospheric pressure (dummy)
1720                           &                 fsdept(ji,jj,jk),             & !      depth
1721                           &                 gphit(ji,jj),                 & !      latitude
1722                           &                 1.0,                          & !      gas transfer (dummy)
1723                           &                 1.0,                          & !      atmospheric xCO2 (dummy)
1724                           &                 1,                            & !      number of points
1725                           &                 PH_OUT,                       & ! OUT: pH
1726                           &                 dummy_out(1),  dummy_out(2),  & !      stuff we don't care about
1727                           &                 dummy_out(3),  dummy_out(4),  &
1728                           &                 dummy_out(5),  dummy_out(6),  &
1729                           &                 dummy_out(7),  dummy_out(8),  &
1730                           &                 dummy_out(9),  dummy_out(10), &
1731                           &                 dummy_out(11), dummy_out(12), &
1732                           &                 dummy_out(13), dummy_out(14), &
1733                           &                 dummy_out(15), dummy_out(16), &
1734                           &                 dummy_out(17), dummy_out(18), &
1735                           &                 dummy_out(19) )
1736
1737                        IF ( jx == 1 ) THEN
1738                           ph_bkg = PH_OUT
1739                        ELSE IF ( jx == 2 ) THEN
1740                           ph_dic = PH_OUT
1741                        ELSE IF ( jx == 3 ) THEN
1742                           ph_alk = PH_OUT
1743                        ENDIF
1744                     END DO
1745
1746                     dph_ddic = (ph_dic - ph_bkg) / zsearch
1747                     dph_dalk = (ph_alk - ph_bkg) / zsearch
1748                     weight   = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) )
1749
1750                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
1751                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
1752                     
1753                  ENDIF
1754                 
1755               END DO
1756            END DO
1757         END DO
1758
1759      ENDIF
1760     
1761      IF ( ll_asmiau ) THEN
1762
1763         !--------------------------------------------------------------------
1764         ! Incremental Analysis Updating
1765         !--------------------------------------------------------------------
1766
1767         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1768
1769            it = kt - nit000 + 1
1770            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1771            ! note this is not a tendency so should not be divided by rdt
1772
1773            IF(lwp) THEN
1774               WRITE(numout,*) 
1775               WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
1776                  &  kt,' with IAU weight = ', pwgtiau(it)
1777               WRITE(numout,*) '~~~~~~~~~~~~'
1778            ENDIF
1779
1780            ! Update the biogeochemical variables
1781            ! Add directly to trn and trb, rather than to tra, because tra gets
1782            ! reset to zero at the start of trc_stp, called after this routine
1783            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1784               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1785               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1786                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1787               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1788                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1789            END WHERE
1790
1791            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1792            ! which is called at end of model run
1793
1794         ENDIF
1795
1796      ELSEIF ( ll_asmdin ) THEN 
1797
1798         !--------------------------------------------------------------------
1799         ! Direct Initialization
1800         !--------------------------------------------------------------------
1801         
1802         IF ( kt == nitdin_r ) THEN
1803
1804            neuler = 0                    ! Force Euler forward step
1805
1806            ! Initialize the now fields with the background + increment
1807            ! Background currently is what the model is initialised with
1808            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', &
1809               &           ' Background state is taken from model rather than background file' )
1810            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1811               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1812               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1813                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1)
1814               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1815            END WHERE
1816 
1817            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1818            ! which is called at end of model run
1819         ENDIF
1820         !
1821      ENDIF
1822#endif     
1823      !
1824   END SUBROUTINE ph_asm_inc
1825
1826   !!===========================================================================
1827   !!===========================================================================
1828   !!===========================================================================
1829
1830   SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1831      !!----------------------------------------------------------------------
1832      !!                    ***  ROUTINE dyn_asm_inc  ***
1833      !!         
1834      !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
1835      !!
1836      !! ** Method  : Direct initialization or Incremental Analysis Updating.
1837      !!
1838      !! ** Action  :
1839      !!----------------------------------------------------------------------
1840      INTEGER,  INTENT(IN) :: kt        ! Current time step
1841      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1842      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1843      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1844      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1845      !
1846      INTEGER  :: jk              ! Loop counter
1847      INTEGER  :: it              ! Index
1848      REAL(wp) :: zincwgt         ! IAU weight for current time step
1849      REAL(wp) :: zincper         ! IAU interval in seconds
1850      !!----------------------------------------------------------------------
1851
1852      IF ( kt <= nit000 ) THEN
1853
1854         !----------------------------------------------------------------------
1855         ! Remove any other balancing increments
1856         !----------------------------------------------------------------------
1857
1858         IF ( ln_pno3inc ) THEN
1859#if defined key_hadocc || defined key_medusa
1860#if defined key_hadocc
1861            it = jp_had_nut
1862#elif defined key_medusa
1863            it = jpdin
1864#endif
1865            IF ( ln_phytobal ) THEN
1866               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1867            ENDIF
1868            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1869               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1870            ENDIF
1871            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1872               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1873            ENDIF
1874            IF ( ln_pphinc ) THEN
1875               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1876            ENDIF
1877#else
1878            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1879#endif
1880         ENDIF
1881
1882         IF ( ln_psi4inc ) THEN
1883#if defined key_medusa
1884            it = jpsil
1885            IF ( ln_phytobal ) THEN
1886               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1887            ENDIF
1888            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1889               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1890            ENDIF
1891            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1892               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1893            ENDIF
1894            IF ( ln_pphinc ) THEN
1895               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1896            ENDIF
1897#else
1898            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1899#endif
1900         ENDIF
1901
1902         IF ( ln_pdicinc ) THEN
1903#if defined key_hadocc || defined key_medusa
1904#if defined key_hadocc
1905            it = jp_had_dic
1906#elif defined key_medusa
1907            it = jpdic
1908#endif
1909            IF ( ln_phytobal ) THEN
1910               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1911            ENDIF
1912            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1913               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1914            ENDIF
1915            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1916               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1917            ENDIF
1918            IF ( ln_pphinc ) THEN
1919               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1920            ENDIF
1921#else
1922            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1923#endif
1924         ENDIF
1925
1926         IF ( ln_palkinc ) THEN
1927#if defined key_hadocc || defined key_medusa
1928#if defined key_hadocc
1929            it = jp_had_alk
1930#elif defined key_medusa
1931            it = jpalk
1932#endif
1933            IF ( ln_phytobal ) THEN
1934               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1935            ENDIF
1936            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1937               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1938            ENDIF
1939            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1940               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1941            ENDIF
1942            IF ( ln_pphinc ) THEN
1943               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1944            ENDIF
1945#else
1946            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1947#endif
1948         ENDIF
1949
1950         IF ( ln_po2inc ) THEN
1951#if defined key_medusa
1952            it = jpoxy
1953            IF ( ln_phytobal ) THEN
1954               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1955            ENDIF
1956            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1957               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1958            ENDIF
1959            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1960               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1961            ENDIF
1962            IF ( ln_pphinc ) THEN
1963               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1964            ENDIF
1965#else
1966            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1967#endif
1968         ENDIF
1969
1970      ENDIF
1971
1972      IF ( ll_asmiau ) THEN
1973
1974         !--------------------------------------------------------------------
1975         ! Incremental Analysis Updating
1976         !--------------------------------------------------------------------
1977
1978         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1979
1980            it = kt - nit000 + 1
1981            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1982            ! note this is not a tendency so should not be divided by rdt
1983
1984            IF(lwp) THEN
1985               WRITE(numout,*) 
1986               WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
1987                  &  kt,' with IAU weight = ', pwgtiau(it)
1988               WRITE(numout,*) '~~~~~~~~~~~~'
1989            ENDIF
1990
1991            ! Update the 3D BGC variables
1992            ! Add directly to trn and trb, rather than to tra, because tra gets
1993            ! reset to zero at the start of trc_stp, called after this routine
1994            ! Don't apply increments if they'll take concentrations negative
1995
1996            IF ( ln_pno3inc ) THEN
1997#if defined key_hadocc
1998               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
1999                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2000                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2001                  trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2002               END WHERE
2003#elif defined key_medusa
2004               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2005                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2006                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2007                  trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2008               END WHERE
2009#else
2010               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2011#endif
2012            ENDIF
2013
2014            IF ( ln_psi4inc ) THEN
2015#if defined key_medusa
2016               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2017                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2018                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2019                  trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2020               END WHERE
2021#else
2022               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2023#endif
2024            ENDIF
2025
2026            IF ( ln_pdicinc ) THEN
2027#if defined key_hadocc
2028               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2029                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2030                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2031                  trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2032               END WHERE
2033#elif defined key_medusa
2034               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2035                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2036                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2037                  trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2038               END WHERE
2039#else
2040               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2041#endif
2042            ENDIF
2043
2044            IF ( ln_palkinc ) THEN
2045#if defined key_hadocc
2046               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2047                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2048                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2049                  trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2050               END WHERE
2051#elif defined key_medusa
2052               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2053                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2054                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2055                  trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2056               END WHERE
2057#else
2058               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2059#endif
2060            ENDIF
2061
2062            IF ( ln_po2inc ) THEN
2063#if defined key_medusa
2064               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2065                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
2066                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2067                  trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2068               END WHERE
2069#else
2070               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2071#endif
2072            ENDIF
2073           
2074            IF ( kt == nitiaufin_r ) THEN
2075               IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2076               IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2077               IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2078               IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2079               IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
2080            ENDIF
2081
2082         ENDIF
2083
2084      ELSEIF ( ll_asmdin ) THEN 
2085
2086         !--------------------------------------------------------------------
2087         ! Direct Initialization
2088         !--------------------------------------------------------------------
2089         
2090         IF ( kt == nitdin_r ) THEN
2091
2092            neuler = 0                    ! Force Euler forward step
2093
2094            ! Initialize the now fields with the background + increment
2095            ! Background currently is what the model is initialised with
2096#if defined key_hadocc
2097            CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
2098               &           ' Background state is taken from model rather than background file' )
2099#elif defined key_medusa
2100            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
2101               &           ' Background state is taken from model rather than background file' )
2102#endif
2103
2104            IF ( ln_pno3inc ) THEN
2105#if defined key_hadocc
2106               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2107                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
2108                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
2109                  trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
2110               END WHERE
2111#elif defined key_medusa
2112               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2113                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
2114                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
2115                  trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
2116               END WHERE
2117#else
2118               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2119#endif
2120            ENDIF
2121
2122            IF ( ln_psi4inc ) THEN
2123#if defined key_medusa
2124               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2125                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
2126                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
2127                  trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
2128               END WHERE
2129#else
2130               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2131#endif
2132            ENDIF
2133
2134            IF ( ln_pdicinc ) THEN
2135#if defined key_hadocc
2136               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2137                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
2138                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
2139                  trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
2140               END WHERE
2141#elif defined key_medusa
2142               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2143                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
2144                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
2145                  trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
2146               END WHERE
2147#else
2148               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2149#endif
2150            ENDIF
2151
2152            IF ( ln_palkinc ) THEN
2153#if defined key_hadocc
2154               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2155                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
2156                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
2157                  trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
2158               END WHERE
2159#elif defined key_medusa
2160               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2161                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
2162                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
2163                  trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
2164               END WHERE
2165#else
2166               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2167#endif
2168            ENDIF
2169
2170            IF ( ln_po2inc ) THEN
2171#if defined key_medusa
2172               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2173                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
2174                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
2175                  trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
2176               END WHERE
2177#else
2178               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2179#endif
2180            ENDIF
2181 
2182            IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2183            IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2184            IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2185            IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2186            IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
2187         ENDIF
2188         !
2189      ENDIF
2190      !
2191   END SUBROUTINE bgc3d_asm_inc
2192
2193   !!===========================================================================
2194
2195END MODULE asmbgc
Note: See TracBrowser for help on using the repository browser.