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

Last change on this file since 11101 was 11101, checked in by frrh, 18 months ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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