source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_phytobal3d/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 13315

Last change on this file since 13315 was 13315, checked in by dford, 4 months ago

Fix for building physics-only.

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