source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

File size: 102.4 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=*),             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=*),                 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      IF(lwp .AND. lflush) CALL flush(numout)
771                                 
772   END SUBROUTINE asm_bgc_bal_wri
773
774   !!===========================================================================
775   !!===========================================================================
776   !!===========================================================================
777
778   SUBROUTINE asm_bgc_bkg_wri( kt, knum )
779      !!------------------------------------------------------------------------
780      !!                    ***  ROUTINE asm_bgc_bkg_wri  ***
781      !!
782      !! ** Purpose :   write out bgc background
783      !!
784      !! ** Method  :   write out bgc background
785      !!
786      !! ** Action  :   write out bgc background
787      !!
788      !! References :   asm_bkg_wri
789      !!------------------------------------------------------------------------
790      !!
791      INTEGER, INTENT(in   ) :: kt     ! Current time-step
792      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
793      !!
794      !!------------------------------------------------------------------------
795
796#if defined key_hadocc
797      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg             )
798      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg             )
799      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg              )
800      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max               )
801      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) )
802      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) )
803      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) )
804      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) )
805      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) )
806      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) )
807      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,:)     )
808      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,:)         )
809#elif defined key_medusa
810      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg        )
811      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        )
812      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         )
813      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          )
814      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) )
815      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd'  , trn(:,:,:,jpchd) )
816      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn'  , trn(:,:,:,jpphn) )
817      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd'  , trn(:,:,:,jpphd) )
818      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds'  , trn(:,:,:,jppds) )
819      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi'  , trn(:,:,:,jpzmi) )
820      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme'  , trn(:,:,:,jpzme) )
821      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din'  , trn(:,:,:,jpdin) )
822      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil'  , trn(:,:,:,jpsil) )
823      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer'  , trn(:,:,:,jpfer) )
824      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det'  , trn(:,:,:,jpdet) )
825      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc'  , trn(:,:,:,jpdtc) )
826      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic'  , trn(:,:,:,jpdic) )
827      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk'  , trn(:,:,:,jpalk) )
828      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
829#endif
830
831   END SUBROUTINE asm_bgc_bkg_wri
832
833   !!===========================================================================
834   !!===========================================================================
835   !!===========================================================================
836
837   SUBROUTINE asm_bgc_unlog_2d( pbkg, pinc_log, pinc_nonlog )
838      !!------------------------------------------------------------------------
839      !!                    ***  ROUTINE asm_bgc_init_incs  ***
840      !!
841      !! ** Purpose :   convert log increments to non-log
842      !!
843      !! ** Method  :   need to account for model background,
844      !!                cannot simply do 10^log_inc. Need to:
845      !!                1) Add log_inc to log10(background) to get log10(analysis)
846      !!                2) Take 10^log10(analysis) to get analysis
847      !!                3) Subtract background from analysis to get non-log incs
848      !!
849      !! ** Action  :   return non-log increments
850      !!
851      !! References :   
852      !!------------------------------------------------------------------------
853      !!
854      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pbkg        ! Background
855      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) :: pinc_log    ! Log incs
856      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) :: pinc_nonlog ! Non-log incs
857      !
858      INTEGER                                     :: ji, jj      ! Loop counters
859      !!
860      !!------------------------------------------------------------------------
861
862      DO jj = 1, jpj
863         DO ji = 1, jpi
864            IF ( pbkg(ji,jj) > 0.0 ) THEN
865               pinc_nonlog(ji,jj) = 10**( LOG10( pbkg(ji,jj) ) + &
866                  &                       pinc_log(ji,jj) )      &
867                  &                 - pbkg(ji,jj)
868            ELSE
869               pinc_nonlog(ji,jj) = 0.0
870            ENDIF
871         END DO
872      END DO
873
874   END SUBROUTINE asm_bgc_unlog_2d
875
876   !!===========================================================================
877   !!===========================================================================
878   !!===========================================================================
879
880   SUBROUTINE phyto2d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
881      !!------------------------------------------------------------------------
882      !!                    ***  ROUTINE phyto2d_asm_inc  ***
883      !!         
884      !! ** Purpose : Apply the chlorophyll assimilation increments.
885      !!
886      !! ** Method  : Calculate increments to state variables using nitrogen
887      !!              balancing.
888      !!              Direct initialization or Incremental Analysis Updating.
889      !!
890      !! ** Action  :
891      !!------------------------------------------------------------------------
892      INTEGER,  INTENT(IN) :: kt        ! Current time step
893      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
894      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
895      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
896      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
897      !
898      INTEGER                      :: it          ! Index
899      REAL(wp)                     :: zincwgt     ! IAU weight for current time step
900      REAL(wp)                     :: zincper     ! IAU interval in seconds
901      REAL(wp), DIMENSION(jpi,jpj) :: zmld        ! Mixed layer depth
902      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chltot ! Local chltot incs
903      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chltot ! Local chltot bkg
904      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phytot ! Local phytot incs
905      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phytot ! Local phytot bkg
906#if defined key_medusa
907      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chldia ! Local chldia incs
908      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chldia ! Local chldia bkg
909      REAL(wp), DIMENSION(jpi,jpj) :: zinc_chlnon ! Local chlnon incs
910      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_chlnon ! Local chlnon bkg
911      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phydia ! Local phydia incs
912      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phydia ! Local phydia bkg
913      REAL(wp), DIMENSION(jpi,jpj) :: zinc_phynon ! Local phynon incs
914      REAL(wp), DIMENSION(jpi,jpj) :: zbkg_phynon ! Local phynon bkg
915#endif
916      !!------------------------------------------------------------------------
917     
918      IF ( kt <= nit000 ) THEN
919
920         ! Un-log any log increments for passing to balancing routines
921         ! Remember that two sets of non-log increments should not be
922         ! expected to be in the same ratio as their log equivalents
923         
924         ! Total chlorophyll
925         IF ( ln_slchltotinc ) THEN
926#if defined key_medusa
927            zbkg_chltot(:,:) = tracer_bkg(:,:,1,jpchn) + tracer_bkg(:,:,1,jpchd)
928#elif defined key_hadocc
929            zbkg_chltot(:,:) = chl_bkg(:,:,1)
930#endif
931            CALL asm_bgc_unlog_2d( zbkg_chltot, slchltot_bkginc, zinc_chltot )
932         ELSE IF ( ln_schltotinc ) THEN
933            zinc_chltot(:,:) = schltot_bkginc(:,:)
934         ELSE
935            zinc_chltot(:,:) = 0.0
936         ENDIF
937
938#if defined key_medusa
939         ! Diatom chlorophyll
940         IF ( ln_slchldiainc ) THEN
941            zbkg_chldia(:,:) = tracer_bkg(:,:,1,jpchd)
942            CALL asm_bgc_unlog_2d( zbkg_chldia, slchldia_bkginc, zinc_chldia )
943         ELSE
944            zinc_chldia(:,:) = 0.0
945         ENDIF
946#endif
947
948#if defined key_medusa
949         ! Non-diatom chlorophyll
950         IF ( ln_slchlnoninc ) THEN
951            zbkg_chlnon(:,:) = tracer_bkg(:,:,1,jpchn)
952            CALL asm_bgc_unlog_2d( zbkg_chlnon, slchlnon_bkginc, zinc_chlnon )
953         ELSE
954            zinc_chlnon(:,:) = 0.0
955         ENDIF
956#endif
957
958         ! Total phytoplankton carbon
959         IF ( ln_slphytotinc ) THEN
960#if defined key_medusa
961            zbkg_phytot(:,:) = (trn(:,:,1,jpphn) * xthetapn) + (trn(:,:,1,jpphd) * xthetapd)
962#elif defined key_hadocc
963            zbkg_phytot(:,:) = trn(:,:,1,jp_had_phy) * c2n_p
964#endif
965            CALL asm_bgc_unlog_2d( zbkg_phytot, slphytot_bkginc, zinc_phytot )
966         ELSE
967            zinc_phytot(:,:) = 0.0
968         ENDIF
969
970#if defined key_medusa
971         ! Diatom phytoplankton carbon
972         IF ( ln_slphydiainc ) THEN
973            zbkg_phydia(:,:) = trn(:,:,1,jpphd) * xthetapd
974            CALL asm_bgc_unlog_2d( zbkg_phydia, slphydia_bkginc, zinc_phydia )
975         ELSE
976            zinc_phydia(:,:) = 0.0
977         ENDIF
978#endif
979
980#if defined key_medusa
981         ! Non-diatom phytoplankton carbon
982         IF ( ln_slphynoninc ) THEN
983            zbkg_phynon(:,:) = trn(:,:,1,jpphn) * xthetapn
984            CALL asm_bgc_unlog_2d( zbkg_phynon, slphynon_bkginc, zinc_phynon )
985         ELSE
986            zinc_phynon(:,:) = 0.0
987         ENDIF
988#endif
989
990         ! Select mixed layer
991         IF ( ll_asmdin ) THEN
992#if defined key_top && ( defined key_hadocc || defined key_medusa )
993            CALL ctl_warn( ' Doing direct initialisation with ocean colour assimilation', &
994               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
995            zmld(:,:) = mld_max_bkg(:,:)
996#else
997            CALL ctl_stop( ' Should not be doing biogeochemical assimilation without key_top' )
998#endif
999         ELSE
1000            SELECT CASE( mld_choice_bgc )
1001            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
1002               zmld(:,:) = hmld(:,:)
1003            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1004               zmld(:,:) = hmlp(:,:)
1005            CASE ( 3 )                   ! Kara MLD [Interpolated]
1006#if defined key_karaml
1007               IF ( ln_kara ) THEN
1008                  zmld(:,:) = hmld_kara(:,:)
1009               ELSE
1010                  CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1011                     &           ' but ln_kara=.false.' )
1012               ENDIF
1013#else
1014               CALL ctl_stop( ' Kara mixed layer requested for phyto2d assimilation,', &
1015                  &           ' but is not defined' )
1016#endif
1017            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
1018               !zmld(:,:) = hmld_tref(:,:)
1019               CALL ctl_stop( ' hmld_tref mixed layer requested for phyto2d assimilation,', &
1020                  &           ' but is not available in this version' )
1021            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1022               zmld(:,:) = hmlpt(:,:)
1023            END SELECT
1024         ENDIF
1025
1026         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
1027
1028#if defined key_medusa
1029         CALL asm_phyto2d_bal_medusa( (ln_slchltotinc .OR. ln_schltotinc), &
1030            &                         zinc_chltot,                         &
1031            &                         ln_slchldiainc,                      &
1032            &                         zinc_chldia,                         &
1033            &                         ln_slchlnoninc,                      &
1034            &                         zinc_chlnon,                         &
1035            &                         ln_slphytotinc,                      &
1036            &                         zinc_phytot,                         &
1037            &                         ln_slphydiainc,                      &
1038            &                         zinc_phydia,                         &
1039            &                         ln_slphynoninc,                      &
1040            &                         zinc_phynon,                         &
1041            &                         zincper,                             &
1042            &                         rn_maxchlinc, ln_phytobal, zmld,     &
1043            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
1044            &                         phyt_avg_bkg, mld_max_bkg,           &
1045            &                         tracer_bkg, phyto2d_balinc )
1046#elif defined key_hadocc
1047         CALL asm_phyto2d_bal_hadocc( (ln_slchltotinc .OR. ln_schltotinc), &
1048            &                         zinc_chltot,                         &
1049            &                         ln_slphytotinc,                      &
1050            &                         zinc_phytot,                         &
1051            &                         zincper,                             &
1052            &                         rn_maxchlinc, ln_phytobal, zmld,     &
1053            &                         pgrow_avg_bkg, ploss_avg_bkg,        &
1054            &                         phyt_avg_bkg, mld_max_bkg,           &
1055            &                         cchl_p_bkg(:,:,1),                   &
1056            &                         tracer_bkg, phyto2d_balinc )
1057#else
1058         CALL ctl_stop( 'Attempting to assimilate phyto2d data, ', &
1059            &           'but not defined a biogeochemical model' )
1060#endif
1061
1062      ENDIF
1063
1064      IF ( ll_asmiau ) THEN
1065
1066         !--------------------------------------------------------------------
1067         ! Incremental Analysis Updating
1068         !--------------------------------------------------------------------
1069
1070         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1071
1072            it = kt - nit000 + 1
1073            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1074            ! note this is not a tendency so should not be divided by rdt
1075
1076            IF(lwp) THEN
1077               WRITE(numout,*) 
1078               WRITE(numout,*) 'phyto2d_asm_inc : phyto2d IAU at time step = ', &
1079                  &  kt,' with IAU weight = ', pwgtiau(it)
1080               WRITE(numout,*) '~~~~~~~~~~~~'
1081            ENDIF
1082
1083            ! Update the biogeochemical variables
1084            ! Add directly to trn and trb, rather than to tra, because tra gets
1085            ! reset to zero at the start of trc_stp, called after this routine
1086#if defined key_medusa
1087            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1088               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1089               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1090                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1091               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1092                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1093            END WHERE
1094#elif defined key_hadocc
1095            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1096               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1097               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1098                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1099               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1100                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1101            END WHERE
1102#endif
1103
1104            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1105            ! which is called at end of model run
1106         ENDIF
1107
1108      ELSEIF ( ll_asmdin ) THEN 
1109
1110         !--------------------------------------------------------------------
1111         ! Direct Initialization
1112         !--------------------------------------------------------------------
1113         
1114         IF ( kt == nitdin_r ) THEN
1115
1116            neuler = 0                    ! Force Euler forward step
1117
1118            ! Initialize the now fields with the background + increment
1119            ! Background currently is what the model is initialised with
1120            CALL ctl_warn( ' Doing direct initialisation with phyto2d assimilation', &
1121               &           ' Background state is taken from model rather than background file' )
1122#if defined key_medusa
1123            WHERE( phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1124               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto2d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1125               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1126                  &                         phyto2d_balinc(:,:,:,jp_msa0:jp_msa1)
1127               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1128            END WHERE
1129#elif defined key_hadocc
1130            WHERE( phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1131               &   trn(:,:,:,jp_had0:jp_had1) + phyto2d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1132               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1133                  &                         phyto2d_balinc(:,:,:,jp_had0:jp_had1)
1134               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1135            END WHERE
1136#endif
1137 
1138            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1139            ! which is called at end of model run
1140         ENDIF
1141         !
1142      ENDIF
1143      !
1144      IF(lwp .AND. lflush) CALL flush(numout)
1145   END SUBROUTINE phyto2d_asm_inc
1146
1147   !!===========================================================================
1148   !!===========================================================================
1149   !!===========================================================================
1150
1151   SUBROUTINE phyto3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1152      !!------------------------------------------------------------------------
1153      !!                    ***  ROUTINE phyto3d_asm_inc  ***
1154      !!         
1155      !! ** Purpose : Apply the profile chlorophyll assimilation increments.
1156      !!
1157      !! ** Method  : Calculate increments to state variables.
1158      !!              Direct initialization or Incremental Analysis Updating.
1159      !!
1160      !! ** Action  :
1161      !!------------------------------------------------------------------------
1162      INTEGER,  INTENT(IN) :: kt        ! Current time step
1163      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1164      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1165      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1166      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1167      !
1168      INTEGER                          :: ji, jj, jk   ! Loop counters
1169      INTEGER                          :: it           ! Index
1170      REAL(wp)                         :: zincwgt      ! IAU weight for timestep
1171      REAL(wp)                         :: zfrac_chn    ! Fraction of jpchn
1172      REAL(wp)                         :: zfrac_chd    ! Fraction of jpchd
1173      REAL(wp)                         :: zrat_phn_chn ! jpphn:jpchn ratio
1174      REAL(wp)                         :: zrat_phd_chd ! jpphd:jpchd ratio
1175      REAL(wp)                         :: zrat_pds_chd ! jppds:jpchd ratio
1176      REAL(wp), DIMENSION(jpi,jpj,jpk) :: chl_inc      ! Chlorophyll increments
1177      REAL(wp), DIMENSION(jpi,jpj,jpk) :: bkg_chl      ! Chlorophyll background
1178      !!------------------------------------------------------------------------
1179
1180      IF ( kt <= nit000 ) THEN
1181
1182         IF ( ln_plchltotinc ) THEN
1183            ! Convert log10(chlorophyll) increment back to a chlorophyll increment
1184            ! In order to transform logchl incs to chl incs, need to account for model
1185            ! background, cannot simply do 10^logchl_bkginc. Need to:
1186            ! 1) Add logchl inc to log10(background) to get log10(analysis)
1187            ! 2) Take 10^log10(analysis) to get analysis
1188            ! 3) Subtract background from analysis to get chl incs
1189            ! If rn_maxchlinc > 0 then cap total absolute chlorophyll increment at that value
1190#if defined key_medusa
1191            bkg_chl(:,:,:) = tracer_bkg(:,:,:,jpchn) + tracer_bkg(:,:,:,jpchd)
1192#elif defined key_hadocc
1193            bkg_chl(:,:,:) = chl_bkg(:,:,:)
1194#endif
1195            DO jk = 1, jpk
1196               DO jj = 1, jpj
1197                  DO ji = 1, jpi
1198                     IF ( bkg_chl(ji,jj,jk) > 0.0 ) THEN
1199                        chl_inc(ji,jj,jk) = 10**( LOG10( bkg_chl(ji,jj,jk) ) + plchltot_bkginc(ji,jj,jk) ) - bkg_chl(ji,jj,jk)
1200                        IF ( rn_maxchlinc > 0.0 ) THEN
1201                           chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( chl_inc(ji,jj,jk), rn_maxchlinc ) )
1202                        ENDIF
1203                     ELSE
1204                        chl_inc(ji,jj,jk) = 0.0
1205                     ENDIF
1206                  END DO
1207               END DO
1208            END DO
1209         ELSE IF ( ln_pchltotinc ) THEN
1210            DO jk = 1, jpk
1211               DO jj = 1, jpj
1212                  DO ji = 1, jpi
1213                     IF ( rn_maxchlinc > 0.0 ) THEN
1214                        chl_inc(ji,jj,jk) = MAX( -1.0 * rn_maxchlinc, MIN( pchltot_bkginc(ji,jj,jk), rn_maxchlinc ) )
1215                     ELSE
1216                        chl_inc(ji,jj,jk) = pchltot_bkginc(ji,jj,jk)
1217                     ENDIF
1218                  END DO
1219               END DO
1220            END DO
1221         ENDIF
1222
1223#if defined key_medusa
1224         ! Loop over each grid point partioning the increments based on existing ratios
1225         DO jk = 1, jpk
1226            DO jj = 1, jpj
1227               DO ji = 1, jpi
1228                  IF ( ( tracer_bkg(ji,jj,jk,jpchn) > 0.0 ) .AND. ( tracer_bkg(ji,jj,jk,jpchd) > 0.0 ) ) THEN
1229                     zfrac_chn = tracer_bkg(ji,jj,jk,jpchn) / (tracer_bkg(ji,jj,jk,jpchn) + tracer_bkg(ji,jj,jk,jpchd))
1230                     zfrac_chd = 1.0 - zfrac_chn
1231                     phyto3d_balinc(ji,jj,jk,jpchn) = chl_inc(ji,jj,jk) * zfrac_chn
1232                     phyto3d_balinc(ji,jj,jk,jpchd) = chl_inc(ji,jj,jk) * zfrac_chd
1233                     zrat_phn_chn = tracer_bkg(ji,jj,jk,jpphn) / tracer_bkg(ji,jj,jk,jpchn)
1234                     zrat_phd_chd = tracer_bkg(ji,jj,jk,jpphd) / tracer_bkg(ji,jj,jk,jpchd)
1235                     zrat_pds_chd = tracer_bkg(ji,jj,jk,jppds) / tracer_bkg(ji,jj,jk,jpchd)
1236                     phyto3d_balinc(ji,jj,jk,jpphn) = phyto3d_balinc(ji,jj,jk,jpchn) * zrat_phn_chn
1237                     phyto3d_balinc(ji,jj,jk,jpphd) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_phd_chd
1238                     phyto3d_balinc(ji,jj,jk,jppds) = phyto3d_balinc(ji,jj,jk,jpchd) * zrat_pds_chd
1239                  ENDIF
1240               END DO
1241            END DO
1242         END DO
1243#elif defined key_hadocc
1244         phyto3d_balinc(:,:,:,jp_had_phy) = ( cchl_p_bkg(:,:,:) / (mw_carbon * c2n_p) ) * chl_inc(:,:,:)
1245#else
1246         CALL ctl_stop( 'Attempting to assimilate p(l)chltot, ', &
1247            &           'but not defined a biogeochemical model' )
1248#endif
1249
1250      ENDIF
1251
1252      IF ( ll_asmiau ) THEN
1253
1254         !--------------------------------------------------------------------
1255         ! Incremental Analysis Updating
1256         !--------------------------------------------------------------------
1257
1258         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1259
1260            it = kt - nit000 + 1
1261            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1262            ! note this is not a tendency so should not be divided by rdt
1263
1264            IF(lwp) THEN
1265               WRITE(numout,*) 
1266               WRITE(numout,*) 'phyto3d_asm_inc : phyto3d IAU at time step = ', &
1267                  &  kt,' with IAU weight = ', pwgtiau(it)
1268               WRITE(numout,*) '~~~~~~~~~~~~'
1269            ENDIF
1270
1271            ! Update the biogeochemical variables
1272            ! Add directly to trn and trb, rather than to tra, because tra gets
1273            ! reset to zero at the start of trc_stp, called after this routine
1274#if defined key_medusa
1275            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1276               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1277               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1278                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1279               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1280                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1281            END WHERE
1282#elif defined key_hadocc
1283            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1284               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1285               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1286                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1287               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1288                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1289            END WHERE
1290#endif
1291
1292            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1293            ! which is called at end of model run
1294         ENDIF
1295
1296      ELSEIF ( ll_asmdin ) THEN 
1297
1298         !--------------------------------------------------------------------
1299         ! Direct Initialization
1300         !--------------------------------------------------------------------
1301         
1302         IF ( kt == nitdin_r ) THEN
1303
1304            neuler = 0                    ! Force Euler forward step
1305
1306            ! Initialize the now fields with the background + increment
1307            ! Background currently is what the model is initialised with
1308            CALL ctl_warn( ' Doing direct initialisation with phyto3d assimilation', &
1309               &           ' Background state is taken from model rather than background file' )
1310#if defined key_medusa
1311            WHERE( phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1312               &   trn(:,:,:,jp_msa0:jp_msa1) + phyto3d_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1313               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1314                  &                         phyto3d_balinc(:,:,:,jp_msa0:jp_msa1)
1315               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1316            END WHERE
1317#elif defined key_hadocc
1318            WHERE( phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1319               &   trn(:,:,:,jp_had0:jp_had1) + phyto3d_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1320               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1321                  &                         phyto3d_balinc(:,:,:,jp_had0:jp_had1)
1322               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1323            END WHERE
1324#endif
1325 
1326            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1327            ! which is called at end of model run
1328         ENDIF
1329         !
1330      ENDIF
1331      !
1332      IF(lwp .AND. lflush) CALL flush(numout)
1333   END SUBROUTINE phyto3d_asm_inc
1334
1335   !!===========================================================================
1336   !!===========================================================================
1337   !!===========================================================================
1338
1339   SUBROUTINE pco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1340      &                     ll_trainc, pt_bkginc, ps_bkginc )
1341      !!------------------------------------------------------------------------
1342      !!                    ***  ROUTINE pco2_asm_inc  ***
1343      !!         
1344      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
1345      !!
1346      !! ** Method  : Calculate increments to state variables using carbon
1347      !!              balancing.
1348      !!              Direct initialization or Incremental Analysis Updating.
1349      !!
1350      !! ** Action  :
1351      !!------------------------------------------------------------------------
1352      INTEGER, INTENT(IN)                   :: kt           ! Current time step
1353      LOGICAL, INTENT(IN)                   :: ll_asmdin    ! Flag for direct initialisation
1354      LOGICAL, INTENT(IN)                   :: ll_asmiau    ! Flag for incremental analysis update
1355      INTEGER, INTENT(IN)                   :: kcycper      ! Dimension of pwgtiau
1356      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1357      LOGICAL, INTENT(IN)                   :: ll_trainc    ! Flag for T/S increments
1358      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1359      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1360      !
1361      INTEGER                               :: ji, jj, jk   ! Loop counters
1362      INTEGER                               :: it           ! Index
1363      INTEGER                               :: jkmax        ! Index of mixed layer depth
1364      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
1365      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
1366      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
1367      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
1368      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
1369      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
1370      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
1371      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
1372      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
1373
1374      ! Coefficients for fCO2 to pCO2 conversion
1375      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
1376      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
1377      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
1378      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
1379      REAL(wp)                              :: zcoef_fco2_5 = 2.0
1380      REAL(wp)                              :: zcoef_fco2_6 = 57.7
1381      REAL(wp)                              :: zcoef_fco2_7 = 0.118
1382      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
1383      !!------------------------------------------------------------------------
1384
1385      IF ( kt <= nit000 ) THEN
1386
1387         !----------------------------------------------------------------------
1388         ! DIC and alkalinity balancing
1389         !----------------------------------------------------------------------
1390
1391         IF ( ln_sfco2inc ) THEN
1392#if defined key_medusa && defined key_roam
1393            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
1394            patm(1) = 1.0
1395            phyd(1) = 0.0
1396            DO jj = 1, jpj
1397               DO ji = 1, jpi
1398                  CALL f2pCO2( sfco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
1399               END DO
1400            END DO
1401#else
1402            ! If assimilating fCO2, then convert to pCO2 using temperature
1403            ! See flux_gas.F90 within HadOCC for details of calculation
1404            pco2_bkginc_temp(:,:) = sfco2_bkginc(:,:) /                                                            &
1405               &                    EXP((zcoef_fco2_1                                                            + &
1406               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
1407               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
1408               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
1409               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
1410               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
1411#endif
1412         ELSE
1413            pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
1414         ENDIF
1415
1416         ! Account for physics assimilation if required
1417         IF ( ll_trainc ) THEN
1418            tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
1419            sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
1420         ELSE
1421            tem_bkg_temp(:,:) = tsn(:,:,1,1)
1422            sal_bkg_temp(:,:) = tsn(:,:,1,2)
1423         ENDIF
1424
1425#if defined key_medusa
1426         ! Account for phytoplankton balancing if required
1427         IF ( ln_phytobal ) THEN
1428            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + phyto2d_balinc(:,:,1,jpdic)
1429            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + phyto2d_balinc(:,:,1,jpalk)
1430         ELSE
1431            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
1432            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
1433         ENDIF
1434
1435         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1436            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1437            &               pco2_balinc(:,:,1,jpdic), pco2_balinc(:,:,1,jpalk) )
1438
1439#elif defined key_hadocc
1440         ! Account for phytoplankton balancing if required
1441         IF ( ln_phytobal ) THEN
1442            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + phyto2d_balinc(:,:,1,jp_had_dic)
1443            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + phyto2d_balinc(:,:,1,jp_had_alk)
1444         ELSE
1445            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
1446            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
1447         ENDIF
1448
1449         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
1450            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
1451            &               pco2_balinc(:,:,1,jp_had_dic), pco2_balinc(:,:,1,jp_had_alk) )
1452
1453#else
1454         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
1455            &           'but not defined a biogeochemical model' )
1456#endif
1457
1458         ! Select mixed layer
1459         IF ( ll_asmdin ) THEN
1460#if defined key_hadocc || defined key_medusa
1461            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1462               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
1463            zmld(:,:) = mld_max_bkg(:,:)
1464#else
1465            CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
1466#endif
1467         ELSE
1468            SELECT CASE( mld_choice_bgc )
1469            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
1470               zmld(:,:) = hmld(:,:)
1471            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
1472               zmld(:,:) = hmlp(:,:)
1473            CASE ( 3 )                   ! Kara MLD [Interpolated]
1474#if defined key_karaml
1475               IF ( ln_kara ) THEN
1476                  zmld(:,:) = hmld_kara(:,:)
1477               ELSE
1478                  CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1479                     &           ' but ln_kara=.false.' )
1480               ENDIF
1481#else
1482               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
1483                  &           ' but is not defined' )
1484#endif
1485            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
1486               !zmld(:,:) = hmld_tref(:,:)
1487               CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
1488                  &           ' but is not available in this version' )
1489            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
1490               zmld(:,:) = hmlpt(:,:)
1491            END SELECT
1492         ENDIF
1493
1494         ! Propagate through mixed layer
1495         DO jj = 1, jpj
1496            DO ji = 1, jpi
1497               !
1498               jkmax = jpk-1
1499               DO jk = jpk-1, 1, -1
1500                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
1501                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
1502                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
1503                     jkmax = jk
1504                  ENDIF
1505               END DO
1506               !
1507#if defined key_top
1508               DO jk = 2, jkmax
1509                  pco2_balinc(ji,jj,jk,:) = pco2_balinc(ji,jj,1,:)
1510               END DO
1511#endif
1512               !
1513            END DO
1514         END DO
1515
1516      ENDIF
1517
1518      IF ( ll_asmiau ) THEN
1519
1520         !--------------------------------------------------------------------
1521         ! Incremental Analysis Updating
1522         !--------------------------------------------------------------------
1523
1524         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1525
1526            it = kt - nit000 + 1
1527            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1528            ! note this is not a tendency so should not be divided by rdt
1529
1530            IF(lwp) THEN
1531               WRITE(numout,*) 
1532               IF ( ln_spco2inc ) THEN
1533                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
1534                     &  kt,' with IAU weight = ', pwgtiau(it)
1535               ELSE IF ( ln_sfco2inc ) THEN
1536                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
1537                     &  kt,' with IAU weight = ', pwgtiau(it)
1538               ENDIF
1539               WRITE(numout,*) '~~~~~~~~~~~~'
1540            ENDIF
1541
1542            ! Update the biogeochemical variables
1543            ! Add directly to trn and trb, rather than to tra, because tra gets
1544            ! reset to zero at the start of trc_stp, called after this routine
1545#if defined key_medusa
1546            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1547               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1548               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1549                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1550               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1551                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1552            END WHERE
1553#elif defined key_hadocc
1554            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1555               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt > 0.0_wp )
1556               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1557                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1558               trb(:,:,:,jp_had0:jp_had1) = trb(:,:,:,jp_had0:jp_had1) + &
1559                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1) * zincwgt
1560            END WHERE
1561#endif
1562
1563            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1564            ! which is called at end of model run
1565
1566         ENDIF
1567
1568      ELSEIF ( ll_asmdin ) THEN 
1569
1570         !--------------------------------------------------------------------
1571         ! Direct Initialization
1572         !--------------------------------------------------------------------
1573         
1574         IF ( kt == nitdin_r ) THEN
1575
1576            neuler = 0                    ! Force Euler forward step
1577
1578            ! Initialize the now fields with the background + increment
1579            ! Background currently is what the model is initialised with
1580            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
1581               &           ' Background state is taken from model rather than background file' )
1582#if defined key_medusa
1583            WHERE( pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1584               &   trn(:,:,:,jp_msa0:jp_msa1) + pco2_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1585               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1586                  &                         pco2_balinc(:,:,:,jp_msa0:jp_msa1)
1587               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1588            END WHERE
1589#elif defined key_hadocc
1590            WHERE( pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp .OR. &
1591               &   trn(:,:,:,jp_had0:jp_had1) + pco2_balinc(:,:,:,jp_had0:jp_had1) > 0.0_wp )
1592               trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
1593                  &                         pco2_balinc(:,:,:,jp_had0:jp_had1)
1594               trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
1595            END WHERE
1596#endif
1597 
1598            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1599            ! which is called at end of model run
1600         ENDIF
1601         !
1602      ENDIF
1603      !
1604      IF(lwp .AND. lflush) CALL flush(numout)
1605   END SUBROUTINE pco2_asm_inc
1606
1607   !!===========================================================================
1608   !!===========================================================================
1609   !!===========================================================================
1610
1611   SUBROUTINE ph_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
1612      &                   ll_trainc, pt_bkginc, ps_bkginc )
1613      !!------------------------------------------------------------------------
1614      !!                    ***  ROUTINE ph_asm_inc  ***
1615      !!         
1616      !! ** Purpose : Apply the pH assimilation increments.
1617      !!
1618      !! ** Method  : Calculate increments to DIC and alkalinity from pH
1619      !!              Use a similar approach to the pCO2 scheme
1620      !!              Direct initialization or Incremental Analysis Updating.
1621      !!
1622      !! ** Action  :
1623      !!------------------------------------------------------------------------
1624      INTEGER,                          INTENT(IN) :: kt        ! Current time step
1625      LOGICAL,                          INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1626      LOGICAL,                          INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1627      INTEGER,                          INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1628      REAL(wp), DIMENSION(kcycper),     INTENT(IN) :: pwgtiau   ! IAU weights
1629      LOGICAL,                          INTENT(IN) :: ll_trainc ! Flag for T/S increments
1630      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
1631      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
1632     
1633      REAL(wp)                         :: zsearch = 10.0 ! Increment to DIC/alk in pH calculation
1634      REAL(wp)                         :: DIC_IN, ALK_IN ! DIC/alk in pH calculation
1635      REAL(wp)                         :: PH_OUT         ! pH from pH calculation
1636      REAL(wp)                         :: ph_bkg         ! pH from pH calculation
1637      REAL(wp)                         :: ph_dic         ! pH from pH calculation
1638      REAL(wp)                         :: ph_alk         ! pH from pH calculation
1639      REAL(wp)                         :: dph_ddic       ! Change in pH wrt DIC
1640      REAL(wp)                         :: dph_dalk       ! Change in pH wrt alk
1641      REAL(wp)                         :: weight         ! Increment weighting
1642      REAL(wp)                         :: zincwgt        ! IAU weight for current time step
1643      REAL(wp), DIMENSION(jpi,jpj,jpk) :: dic_bkg_temp   ! DIC background
1644      REAL(wp), DIMENSION(jpi,jpj,jpk) :: alk_bkg_temp   ! Alkalinity background
1645      REAL(wp), DIMENSION(jpi,jpj,jpk) :: din_bkg_temp   ! DIN background
1646      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sil_bkg_temp   ! Silicate background
1647      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tem_bkg_temp   ! Temperature background
1648      REAL(wp), DIMENSION(jpi,jpj,jpk) :: sal_bkg_temp   ! Salinity background
1649      REAL(wp), DIMENSION(19)          :: dummy_out      ! Dummy variables
1650      INTEGER                          :: ji, jj, jk, jx ! Loop counters
1651      INTEGER                          :: it             ! Index
1652      !!------------------------------------------------------------------------
1653
1654#if ! defined key_medusa
1655      CALL ctl_stop( ' pH balancing only implemented for MEDUSA' )
1656#else
1657
1658      IF ( kt <= nit000 ) THEN
1659
1660         !----------------------------------------------------------------------
1661         ! DIC and alkalinity balancing
1662         !----------------------------------------------------------------------
1663
1664         ! Account for physics assimilation if required
1665         IF ( ll_trainc ) THEN
1666            tem_bkg_temp(:,:,:) = tsn(:,:,:,1) + pt_bkginc(:,:,:)
1667            sal_bkg_temp(:,:,:) = tsn(:,:,:,2) + ps_bkginc(:,:,:)
1668         ELSE
1669            tem_bkg_temp(:,:,:) = tsn(:,:,:,1)
1670            sal_bkg_temp(:,:,:) = tsn(:,:,:,2)
1671         ENDIF
1672
1673         ! Account for phytoplankton balancing if required
1674         IF ( ln_phytobal ) THEN
1675            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic) + phyto2d_balinc(:,:,:,jpdic)
1676            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk) + phyto2d_balinc(:,:,:,jpalk)
1677            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin) + phyto2d_balinc(:,:,:,jpdin)
1678            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil) + phyto2d_balinc(:,:,:,jpsil)
1679         ELSE
1680            dic_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdic)
1681            alk_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpalk)
1682            din_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpdin)
1683            sil_bkg_temp(:,:,:) = tracer_bkg(:,:,:,jpsil)
1684         ENDIF
1685         
1686         ! Account for pCO2 balancing if required
1687         IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1688            dic_bkg_temp(:,:,:) = dic_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpdic)
1689            alk_bkg_temp(:,:,:) = alk_bkg_temp(:,:,:) + pco2_balinc(:,:,:,jpalk)
1690         ENDIF
1691         
1692         ! Loop over grid points and calculate dpH/dDIC and dpH/dAlk
1693         ! This requires three calls to the MOCSY carbonate package
1694         ! One to find pH at background DIC and alkalinity
1695         ! One to find pH when DIC is increased by 10
1696         ! One to find pH when alkalinity is increased by 10
1697         ! Then calculate the gradients
1698
1699         DO jk = 1, jpk
1700            DO jj = 1, jpj
1701               DO ji = 1, jpi
1702
1703                  IF ( pph_bkginc(ji,jj,jk) /= 0.0 ) THEN
1704
1705                     DO jx = 1, 3
1706                        IF ( jx == 1 ) THEN
1707                           DIC_IN = dic_bkg_temp(ji,jj,jk)
1708                           ALK_IN = alk_bkg_temp(ji,jj,jk)
1709                        ELSE IF ( jx == 2 ) THEN
1710                           DIC_IN = dic_bkg_temp(ji,jj,jk) + zsearch
1711                           ALK_IN = alk_bkg_temp(ji,jj,jk)
1712                        ELSE IF ( jx == 3 ) THEN
1713                           DIC_IN = dic_bkg_temp(ji,jj,jk)
1714                           ALK_IN = alk_bkg_temp(ji,jj,jk) + zsearch
1715                        ENDIF
1716
1717                        CALL mocsy_interface(tsn(ji,jj,jk,jp_tem),         & ! IN:  temperature
1718                           &                 tsn(ji,jj,jk,jp_sal),         & !      salinity
1719                           &                 ALK_IN,                       & !      alkalinity
1720                           &                 DIC_IN,                       & !      DIC
1721                           &                 sil_bkg_temp(ji,jj,jk),       & !      silicate
1722                           &                 din_bkg_temp(ji,jj,jk)/16.0,  & !      phosphate (Redfield from DIN)
1723                           &                 1.0,                          & !      atmospheric pressure (dummy)
1724                           &                 fsdept(ji,jj,jk),             & !      depth
1725                           &                 gphit(ji,jj),                 & !      latitude
1726                           &                 1.0,                          & !      gas transfer (dummy)
1727                           &                 1.0,                          & !      atmospheric xCO2 (dummy)
1728                           &                 1,                            & !      number of points
1729                           &                 PH_OUT,                       & ! OUT: pH
1730                           &                 dummy_out(1),  dummy_out(2),  & !      stuff we don't care about
1731                           &                 dummy_out(3),  dummy_out(4),  &
1732                           &                 dummy_out(5),  dummy_out(6),  &
1733                           &                 dummy_out(7),  dummy_out(8),  &
1734                           &                 dummy_out(9),  dummy_out(10), &
1735                           &                 dummy_out(11), dummy_out(12), &
1736                           &                 dummy_out(13), dummy_out(14), &
1737                           &                 dummy_out(15), dummy_out(16), &
1738                           &                 dummy_out(17), dummy_out(18), &
1739                           &                 dummy_out(19) )
1740
1741                        IF ( jx == 1 ) THEN
1742                           ph_bkg = PH_OUT
1743                        ELSE IF ( jx == 2 ) THEN
1744                           ph_dic = PH_OUT
1745                        ELSE IF ( jx == 3 ) THEN
1746                           ph_alk = PH_OUT
1747                        ENDIF
1748                     END DO
1749
1750                     dph_ddic = (ph_dic - ph_bkg) / zsearch
1751                     dph_dalk = (ph_alk - ph_bkg) / zsearch
1752                     weight   = pph_bkginc(ji,jj,jk) / ( (dph_ddic*dph_ddic) + (dph_dalk*dph_dalk) )
1753
1754                     ph_balinc(ji,jj,jk,jpdic) = weight * dph_ddic
1755                     ph_balinc(ji,jj,jk,jpalk) = weight * dph_dalk
1756                     
1757                  ENDIF
1758                 
1759               END DO
1760            END DO
1761         END DO
1762
1763      ENDIF
1764     
1765      IF ( ll_asmiau ) THEN
1766
1767         !--------------------------------------------------------------------
1768         ! Incremental Analysis Updating
1769         !--------------------------------------------------------------------
1770
1771         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1772
1773            it = kt - nit000 + 1
1774            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1775            ! note this is not a tendency so should not be divided by rdt
1776
1777            IF(lwp) THEN
1778               WRITE(numout,*) 
1779               WRITE(numout,*) 'ph_asm_inc : pH IAU at time step = ', &
1780                  &  kt,' with IAU weight = ', pwgtiau(it)
1781               WRITE(numout,*) '~~~~~~~~~~~~'
1782            ENDIF
1783
1784            ! Update the biogeochemical variables
1785            ! Add directly to trn and trb, rather than to tra, because tra gets
1786            ! reset to zero at the start of trc_stp, called after this routine
1787            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1788               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt > 0.0_wp )
1789               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1790                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1791               trb(:,:,:,jp_msa0:jp_msa1) = trb(:,:,:,jp_msa0:jp_msa1) + &
1792                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1) * zincwgt
1793            END WHERE
1794
1795            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1796            ! which is called at end of model run
1797
1798         ENDIF
1799
1800      ELSEIF ( ll_asmdin ) THEN 
1801
1802         !--------------------------------------------------------------------
1803         ! Direct Initialization
1804         !--------------------------------------------------------------------
1805         
1806         IF ( kt == nitdin_r ) THEN
1807
1808            neuler = 0                    ! Force Euler forward step
1809
1810            ! Initialize the now fields with the background + increment
1811            ! Background currently is what the model is initialised with
1812            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pH assimilation', &
1813               &           ' Background state is taken from model rather than background file' )
1814            WHERE( ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp .OR. &
1815               &   trn(:,:,:,jp_msa0:jp_msa1) + ph_balinc(:,:,:,jp_msa0:jp_msa1) > 0.0_wp )
1816               trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
1817                  &                         ph_balinc(:,:,:,jp_msa0:jp_msa1)
1818               trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
1819            END WHERE
1820 
1821            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
1822            ! which is called at end of model run
1823         ENDIF
1824         !
1825      ENDIF
1826#endif     
1827      !
1828      IF(lwp .AND. lflush) CALL flush(numout)
1829
1830   END SUBROUTINE ph_asm_inc
1831
1832   !!===========================================================================
1833   !!===========================================================================
1834   !!===========================================================================
1835
1836   SUBROUTINE bgc3d_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
1837      !!----------------------------------------------------------------------
1838      !!                    ***  ROUTINE dyn_asm_inc  ***
1839      !!         
1840      !! ** Purpose : Apply generic 3D biogeochemistry assimilation increments.
1841      !!
1842      !! ** Method  : Direct initialization or Incremental Analysis Updating.
1843      !!
1844      !! ** Action  :
1845      !!----------------------------------------------------------------------
1846      INTEGER,  INTENT(IN) :: kt        ! Current time step
1847      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
1848      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
1849      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
1850      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
1851      !
1852      INTEGER  :: it              ! Index
1853      REAL(wp) :: zincwgt         ! IAU weight for current time step
1854      !!----------------------------------------------------------------------
1855
1856      IF ( kt <= nit000 ) THEN
1857
1858         !----------------------------------------------------------------------
1859         ! Remove any other balancing increments
1860         !----------------------------------------------------------------------
1861
1862         IF ( ln_pno3inc ) THEN
1863#if defined key_hadocc || defined key_medusa
1864#if defined key_hadocc
1865            it = jp_had_nut
1866#elif defined key_medusa
1867            it = jpdin
1868#endif
1869            IF ( ln_phytobal ) THEN
1870               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1871            ENDIF
1872            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1873               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1874            ENDIF
1875            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1876               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1877            ENDIF
1878            IF ( ln_pphinc ) THEN
1879               pno3_bkginc(:,:,:) = pno3_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1880            ENDIF
1881#else
1882            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1883#endif
1884         ENDIF
1885
1886         IF ( ln_psi4inc ) THEN
1887#if defined key_medusa
1888            it = jpsil
1889            IF ( ln_phytobal ) THEN
1890               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1891            ENDIF
1892            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1893               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1894            ENDIF
1895            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1896               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1897            ENDIF
1898            IF ( ln_pphinc ) THEN
1899               psi4_bkginc(:,:,:) = psi4_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1900            ENDIF
1901#else
1902            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1903#endif
1904         ENDIF
1905
1906         IF ( ln_pdicinc ) THEN
1907#if defined key_hadocc || defined key_medusa
1908#if defined key_hadocc
1909            it = jp_had_dic
1910#elif defined key_medusa
1911            it = jpdic
1912#endif
1913            IF ( ln_phytobal ) THEN
1914               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1915            ENDIF
1916            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1917               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1918            ENDIF
1919            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1920               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1921            ENDIF
1922            IF ( ln_pphinc ) THEN
1923               pdic_bkginc(:,:,:) = pdic_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1924            ENDIF
1925#else
1926            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1927#endif
1928         ENDIF
1929
1930         IF ( ln_palkinc ) THEN
1931#if defined key_hadocc || defined key_medusa
1932#if defined key_hadocc
1933            it = jp_had_alk
1934#elif defined key_medusa
1935            it = jpalk
1936#endif
1937            IF ( ln_phytobal ) THEN
1938               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1939            ENDIF
1940            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1941               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1942            ENDIF
1943            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1944               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1945            ENDIF
1946            IF ( ln_pphinc ) THEN
1947               palk_bkginc(:,:,:) = palk_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1948            ENDIF
1949#else
1950            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1951#endif
1952         ENDIF
1953
1954         IF ( ln_po2inc ) THEN
1955#if defined key_medusa
1956            it = jpoxy
1957            IF ( ln_phytobal ) THEN
1958               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto2d_balinc(:,:,:,it)
1959            ENDIF
1960            IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1961               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - phyto3d_balinc(:,:,:,it)
1962            ENDIF
1963            IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1964               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - pco2_balinc(:,:,:,it)
1965            ENDIF
1966            IF ( ln_pphinc ) THEN
1967               po2_bkginc(:,:,:) = po2_bkginc(:,:,:) - ph_balinc(:,:,:,it)
1968            ENDIF
1969#else
1970            CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
1971#endif
1972         ENDIF
1973
1974      ENDIF
1975
1976      IF ( ll_asmiau ) THEN
1977
1978         !--------------------------------------------------------------------
1979         ! Incremental Analysis Updating
1980         !--------------------------------------------------------------------
1981
1982         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1983
1984            it = kt - nit000 + 1
1985            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
1986            ! note this is not a tendency so should not be divided by rdt
1987
1988            IF(lwp) THEN
1989               WRITE(numout,*) 
1990               WRITE(numout,*) 'bgc3d_asm_inc : 3D BGC IAU at time step = ', &
1991                  &  kt,' with IAU weight = ', pwgtiau(it)
1992               WRITE(numout,*) '~~~~~~~~~~~~'
1993            ENDIF
1994
1995            ! Update the 3D BGC variables
1996            ! Add directly to trn and trb, rather than to tra, because tra gets
1997            ! reset to zero at the start of trc_stp, called after this routine
1998            ! Don't apply increments if they'll take concentrations negative
1999
2000            IF ( ln_pno3inc ) THEN
2001#if defined key_hadocc
2002               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2003                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2004                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2005                  trb(:,:,:,jp_had_nut) = trb(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) * zincwgt
2006               END WHERE
2007#elif defined key_medusa
2008               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2009                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt > 0.0_wp )
2010                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2011                  trb(:,:,:,jpdin) = trb(:,:,:,jpdin) + pno3_bkginc(:,:,:) * zincwgt
2012               END WHERE
2013#else
2014               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2015#endif
2016            ENDIF
2017
2018            IF ( ln_psi4inc ) THEN
2019#if defined key_medusa
2020               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2021                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt > 0.0_wp )
2022                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2023                  trb(:,:,:,jpsil) = trb(:,:,:,jpsil) + psi4_bkginc(:,:,:) * zincwgt
2024               END WHERE
2025#else
2026               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2027#endif
2028            ENDIF
2029
2030            IF ( ln_pdicinc ) THEN
2031#if defined key_hadocc
2032               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2033                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2034                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2035                  trb(:,:,:,jp_had_dic) = trb(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) * zincwgt
2036               END WHERE
2037#elif defined key_medusa
2038               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2039                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt > 0.0_wp )
2040                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2041                  trb(:,:,:,jpdic) = trb(:,:,:,jpdic) + pdic_bkginc(:,:,:) * zincwgt
2042               END WHERE
2043#else
2044               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2045#endif
2046            ENDIF
2047
2048            IF ( ln_palkinc ) THEN
2049#if defined key_hadocc
2050               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2051                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2052                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2053                  trb(:,:,:,jp_had_alk) = trb(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) * zincwgt
2054               END WHERE
2055#elif defined key_medusa
2056               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2057                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt > 0.0_wp )
2058                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2059                  trb(:,:,:,jpalk) = trb(:,:,:,jpalk) + palk_bkginc(:,:,:) * zincwgt
2060               END WHERE
2061#else
2062               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2063#endif
2064            ENDIF
2065
2066            IF ( ln_po2inc ) THEN
2067#if defined key_medusa
2068               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2069                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt > 0.0_wp )
2070                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2071                  trb(:,:,:,jpoxy) = trb(:,:,:,jpoxy) + po2_bkginc(:,:,:) * zincwgt
2072               END WHERE
2073#else
2074               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2075#endif
2076            ENDIF
2077           
2078            IF ( kt == nitiaufin_r ) THEN
2079               IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2080               IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2081               IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2082               IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2083               IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
2084            ENDIF
2085
2086         ENDIF
2087
2088      ELSEIF ( ll_asmdin ) THEN 
2089
2090         !--------------------------------------------------------------------
2091         ! Direct Initialization
2092         !--------------------------------------------------------------------
2093         
2094         IF ( kt == nitdin_r ) THEN
2095
2096            neuler = 0                    ! Force Euler forward step
2097
2098            ! Initialize the now fields with the background + increment
2099            ! Background currently is what the model is initialised with
2100#if defined key_hadocc
2101            CALL ctl_warn( ' Doing direct initialisation of HadOCC with 3D BGC assimilation', &
2102               &           ' Background state is taken from model rather than background file' )
2103#elif defined key_medusa
2104            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with 3D BGC assimilation', &
2105               &           ' Background state is taken from model rather than background file' )
2106#endif
2107
2108            IF ( ln_pno3inc ) THEN
2109#if defined key_hadocc
2110               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2111                  &   trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:) > 0.0_wp )
2112                  trn(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut) + pno3_bkginc(:,:,:)
2113                  trb(:,:,:,jp_had_nut) = trn(:,:,:,jp_had_nut)
2114               END WHERE
2115#elif defined key_medusa
2116               WHERE( pno3_bkginc(:,:,:) > 0.0_wp .OR. &
2117                  &   trn(:,:,:,jpdin) + pno3_bkginc(:,:,:) > 0.0_wp )
2118                  trn(:,:,:,jpdin) = trn(:,:,:,jpdin) + pno3_bkginc(:,:,:)
2119                  trb(:,:,:,jpdin) = trn(:,:,:,jpdin)
2120               END WHERE
2121#else
2122               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2123#endif
2124            ENDIF
2125
2126            IF ( ln_psi4inc ) THEN
2127#if defined key_medusa
2128               WHERE( psi4_bkginc(:,:,:) > 0.0_wp .OR. &
2129                  &   trn(:,:,:,jpsil) + psi4_bkginc(:,:,:) > 0.0_wp )
2130                  trn(:,:,:,jpsil) = trn(:,:,:,jpsil) + psi4_bkginc(:,:,:)
2131                  trb(:,:,:,jpsil) = trn(:,:,:,jpsil)
2132               END WHERE
2133#else
2134               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2135#endif
2136            ENDIF
2137
2138            IF ( ln_pdicinc ) THEN
2139#if defined key_hadocc
2140               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2141                  &   trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:) > 0.0_wp )
2142                  trn(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic) + pdic_bkginc(:,:,:)
2143                  trb(:,:,:,jp_had_dic) = trn(:,:,:,jp_had_dic)
2144               END WHERE
2145#elif defined key_medusa
2146               WHERE( pdic_bkginc(:,:,:) > 0.0_wp .OR. &
2147                  &   trn(:,:,:,jpdic) + pdic_bkginc(:,:,:) > 0.0_wp )
2148                  trn(:,:,:,jpdic) = trn(:,:,:,jpdic) + pdic_bkginc(:,:,:)
2149                  trb(:,:,:,jpdic) = trn(:,:,:,jpdic)
2150               END WHERE
2151#else
2152               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2153#endif
2154            ENDIF
2155
2156            IF ( ln_palkinc ) THEN
2157#if defined key_hadocc
2158               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2159                  &   trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:) > 0.0_wp )
2160                  trn(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk) + palk_bkginc(:,:,:)
2161                  trb(:,:,:,jp_had_alk) = trn(:,:,:,jp_had_alk)
2162               END WHERE
2163#elif defined key_medusa
2164               WHERE( palk_bkginc(:,:,:) > 0.0_wp .OR. &
2165                  &   trn(:,:,:,jpalk) + palk_bkginc(:,:,:) > 0.0_wp )
2166                  trn(:,:,:,jpalk) = trn(:,:,:,jpalk) + palk_bkginc(:,:,:)
2167                  trb(:,:,:,jpalk) = trn(:,:,:,jpalk)
2168               END WHERE
2169#else
2170               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2171#endif
2172            ENDIF
2173
2174            IF ( ln_po2inc ) THEN
2175#if defined key_medusa
2176               WHERE( po2_bkginc(:,:,:) > 0.0_wp .OR. &
2177                  &   trn(:,:,:,jpoxy) + po2_bkginc(:,:,:) > 0.0_wp )
2178                  trn(:,:,:,jpoxy) = trn(:,:,:,jpoxy) + po2_bkginc(:,:,:)
2179                  trb(:,:,:,jpoxy) = trn(:,:,:,jpoxy)
2180               END WHERE
2181#else
2182               CALL ctl_stop ( ' bgc3d_asm_inc: no compatible BGC model defined' )
2183#endif
2184            ENDIF
2185 
2186            IF ( ln_pno3inc ) DEALLOCATE( pno3_bkginc )
2187            IF ( ln_psi4inc ) DEALLOCATE( psi4_bkginc )
2188            IF ( ln_pdicinc ) DEALLOCATE( pdic_bkginc )
2189            IF ( ln_palkinc ) DEALLOCATE( palk_bkginc )
2190            IF ( ln_po2inc  ) DEALLOCATE( po2_bkginc  )
2191         ENDIF
2192         !
2193      ENDIF
2194      !
2195      IF(lwp .AND. lflush) CALL flush(numout)
2196   END SUBROUTINE bgc3d_asm_inc
2197
2198   !!===========================================================================
2199
2200END MODULE asmbgc
Note: See TracBrowser for help on using the repository browser.