source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 13191

Last change on this file since 13191 was 13191, checked in by jwhile, 11 months ago

Updates for 1d runnig

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