source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 9862

Last change on this file since 9862 was 9862, checked in by dford, 2 years ago

Fix to read in DIC and alkalinity background if doing pH assimilation alongside ocean colour assimilation without nitrogen balancing.

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