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

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

Merge in functionality added to GO6 at r10149.

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