source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_v3/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 10057

Last change on this file since 10057 was 10057, checked in by dford, 3 years ago

Fixes for physics-only compilation.

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