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

Last change on this file since 13097 was 13097, checked in by dford, 6 months ago

Allow nitrogen balancing for both 2D and 3D chlorophyll increments.

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