New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asmbgc.F90 in branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asmbgc.F90 @ 9296

Last change on this file since 9296 was 9296, checked in by dford, 6 years ago

Move most of the BGC assimilation code into its own module, to make it more portable across versions/configurations. Partially apply the new naming conventions for BGC obs.

File size: 44.1 KB
Line 
1MODULE asmbgc
2   !!======================================================================
3   !!                       ***  MODULE asmbgc  ***
4   !! Routines and declarations for biogeochemical assimilation
5   !!======================================================================
6   !! History :  3.6  ! 2018-02 (D. Ford)  Adapted from existing modules
7   !!----------------------------------------------------------------------
8   !! 'key_asminc'          : assimilation increment interface
9   !! 'key_top'             : passive tracers
10   !! 'key_hadocc'          : HadOCC model
11   !! 'key_medusa'          : MEDUSA model
12   !! 'key_foam_medusa'     : MEDUSA extras for FOAM OBS and ASM
13   !! 'key_roam'            : MEDUSA extras for carbonate cycle
14   !! 'key_karaml'          : Kara mixed layer depth
15   !!----------------------------------------------------------------------
16   !! asm_bgc_check_options : check biogeochemical assimilation options
17   !! asm_bgc_init_incs     : initialise bgc increments
18   !! asm_bgc_init_bkg      : initialise bgc background
19   !! asm_bgc_bal_wri       : write out balancing increments
20   !! asm_bgc_bkg_wri       : write out bgc background
21   !! slchltot_asm_inc      : apply the slchltot increment
22   !! spco2_asm_inc         : apply the pCO2/fCO2 increment
23   !!----------------------------------------------------------------------
24   USE par_kind,      ONLY: wp             ! kind parameters
25   USE par_oce,       ONLY: jpi, jpj, jpk  ! domain array sizes
26   USE iom                                 ! i/o
27   USE zdfmxl, ONLY :  & 
28!   &  hmld_tref,       &   
29#if defined key_karaml
30   &  hmld_kara,       &
31   &  ln_kara,         &
32#endif   
33   &  hmld,            & 
34   &  hmlp,            &
35   &  hmlpt 
36#if defined key_top
37   USE trc, ONLY: &
38      & trn,      &
39      & trb
40   USE par_trc, ONLY: &
41      & jptra
42#endif
43#if defined key_medusa && defined key_foam_medusa
44   USE asmlogchlbal_medusa, ONLY: &
45      & asm_logchl_bal_medusa
46   USE asmpco2bal, ONLY: &
47      & asm_pco2_bal
48   USE par_medusa
49   USE mocsy_mainmod     ! f2pCO2
50   USE sms_medusa, ONLY: pgrow_avg, &
51      &                  ploss_avg, &
52      &                  phyt_avg,  &
53      &                  mld_max
54#elif defined key_hadocc
55   USE asmlogchlbal_hadocc, ONLY: &
56      & asm_logchl_bal_hadocc
57   USE asmpco2bal, ONLY: &
58      & asm_pco2_bal
59   USE par_hadocc
60   USE trc, ONLY: pgrow_avg, &
61       &     ploss_avg, &
62       &     phyt_avg,  &
63       &     mld_max,   &
64       &     HADOCC_CHL
65   USE had_bgc_const, ONLY: cchl_p
66#endif
67   USE asmpar    ! c_asmbkg, c_asmbal, nitbkg_r, nitdin_r, nitiaustr_r, nitiaufin_r
68   USE oce              ! tsn
69
70   IMPLICIT NONE
71   PRIVATE                   
72
73   PUBLIC asm_bgc_check_options
74   PUBLIC asm_bgc_init_incs
75   PUBLIC asm_bgc_init_bkg
76   PUBLIC asm_bgc_bal_wri
77   PUBLIC asm_bgc_bkg_wri
78   PUBLIC slchltot_asm_inc
79   PUBLIC spco2_asm_inc
80
81   LOGICAL, PUBLIC :: ln_balwri = .FALSE.        !: No output of balancing incs
82   LOGICAL, PUBLIC :: ln_slchltotinc = .FALSE.   !: No slchltot increment
83   LOGICAL, PUBLIC :: ln_slchltotbal = .FALSE.   !: No slchltot balancing
84   LOGICAL, PUBLIC :: ln_spco2inc = .FALSE.      !: No pCO2 increment
85   LOGICAL, PUBLIC :: ln_sfco2inc = .FALSE.      !: No fCO2 increment
86
87   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: slchltot_bkginc !: Increment to background slchltot
88   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: spco2_bkginc   !: Increment to background pCO2
89#if defined key_top
90   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: slchltot_balinc  !: Increment to BGC variables from logchl assim
91   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: spco2_balinc    !: Increment to BGC variables from pCO2 assim
92#endif
93#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
94   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: pgrow_avg_bkg  !: Background phyto growth for logchl balancing
95   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: ploss_avg_bkg  !: Background phyto loss for logchl balancing
96   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: phyt_avg_bkg   !: Background phyto for logchl balancing
97   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: mld_max_bkg    !: Background max MLD for logchl balancing
98   REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: tracer_bkg     !: Background tracer state variables
99#endif
100#if defined key_hadocc
101   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: chl_bkg        !: Background surface chlorophyll
102   REAL(wp), DIMENSION(:,:),     ALLOCATABLE :: cchl_p_bkg     !: Background surface carbon:chlorophyll
103#endif
104   REAL(wp), PUBLIC :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation
105                                      !: <= 0 implies no maximum applied (switch turned off)
106                                      !:  > 0 implies maximum absolute chl increment capped at this value
107
108   INTEGER, PUBLIC :: mld_choice_bgc    = 1   !: choice of mld criteria to use for biogeochemistry assimilation
109                                      !: 1) hmld      - Turbocline/mixing depth                           [W points]
110                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points]
111                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated]
112                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points]
113                                      !: 5) hmlpt     - Density criterion (0.01 kg/m^3 change from 10m)   [T points]
114
115CONTAINS
116
117   SUBROUTINE asm_bgc_check_options
118      !!---------------------------------------------------------------------------
119      !!                    ***  ROUTINE asm_bgc_check_options  ***
120      !!
121      !! ** Purpose :   check validity of biogeochemical assimilation options
122      !!
123      !! ** Method  :   check settings of logicals
124      !!
125      !! ** Action  :   call ctl_stop/ctl_warn if required
126      !!
127      !! References :   asm_inc_init
128      !!---------------------------------------------------------------------------
129
130      IF ( ( ln_balwri ).AND.( .NOT. ln_slchltotinc ).AND.( .NOT. ln_spco2inc ).AND.( .NOT. ln_sfco2inc ) ) THEN
131         CALL ctl_warn( ' Balancing increments are only calculated for logchl and pCO2/fCO2', &
132            &           ' Not assimilating logchl, pCO2 or fCO2, so ln_balwri will be set to .false.')
133         ln_balwri = .FALSE.
134      ENDIF
135
136      IF ( ( ln_slchltotbal ).AND.( .NOT. ln_slchltotinc ) ) THEN
137         CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', &
138            &           ' Not assimilating logchl, so ln_slchltotbal will be set to .false.')
139         ln_slchltotbal = .FALSE.
140      ENDIF
141
142      IF ( ( ln_spco2inc ).AND.( ln_sfco2inc ) ) THEN
143         CALL ctl_stop( ' Can only assimilate pCO2 OR fCO2, not both' )
144      ENDIF
145
146   END SUBROUTINE asm_bgc_check_options
147
148   !!======================================================================
149   !!======================================================================
150   !!======================================================================
151
152   SUBROUTINE asm_bgc_init_incs( knum )
153      !!---------------------------------------------------------------------------
154      !!                    ***  ROUTINE asm_bgc_init_incs  ***
155      !!
156      !! ** Purpose :   initialise bgc increments
157      !!
158      !! ** Method  :   allocate arrays and read increments from file
159      !!
160      !! ** Action  :   allocate arrays and read increments from file
161      !!
162      !! References :   asm_inc_init
163      !!---------------------------------------------------------------------------
164      !!
165      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
166      !!
167      !!---------------------------------------------------------------------------
168
169      IF ( ln_slchltotinc ) THEN
170         ALLOCATE( slchltot_bkginc(jpi,jpj) )
171         slchltot_bkginc(:,:) = 0.0
172#if defined key_top
173         ALLOCATE( slchltot_balinc(jpi,jpj,jpk,jptra) )
174         slchltot_balinc(:,:,:,:) = 0.0
175#endif
176      ENDIF
177      IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
178         ALLOCATE( spco2_bkginc(jpi,jpj) )
179         spco2_bkginc(:,:) = 0.0
180#if defined key_top
181         ALLOCATE( spco2_balinc(jpi,jpj,jpk,jptra) )
182         spco2_balinc(:,:,:,:) = 0.0
183#endif
184      ENDIF
185
186      IF ( ln_slchltotinc ) THEN
187         CALL iom_get( knum, jpdom_autoglo, 'bckinslchltot', slchltot_bkginc(:,:), 1 )
188         ! Apply the masks
189         slchltot_bkginc(:,:) = slchltot_bkginc(:,:) * tmask(:,:,1)
190         ! Set missing increments to 0.0 rather than 1e+20
191         ! to allow for differences in masks
192         WHERE( ABS( slchltot_bkginc(:,:) ) > 1.0e+10 ) slchltot_bkginc(:,:) = 0.0
193      ENDIF
194
195      IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
196         IF ( ln_spco2inc ) THEN
197            CALL iom_get( knum, jpdom_autoglo, 'bckinspco2', spco2_bkginc(:,:), 1 )
198         ELSE IF ( ln_sfco2inc ) THEN
199            CALL iom_get( knum, jpdom_autoglo, 'bckinsfco2', spco2_bkginc(:,:), 1 )
200         ENDIF
201         ! Apply the masks
202         spco2_bkginc(:,:) = spco2_bkginc(:,:) * tmask(:,:,1)
203         ! Set missing increments to 0.0 rather than 1e+20
204         ! to allow for differences in masks
205         WHERE( ABS( spco2_bkginc(:,:) ) > 1.0e+10 ) spco2_bkginc(:,:) = 0.0
206      ENDIF
207
208   END SUBROUTINE asm_bgc_init_incs
209
210   !!======================================================================
211   !!======================================================================
212   !!======================================================================
213
214   SUBROUTINE asm_bgc_init_bkg
215      !!---------------------------------------------------------------------------
216      !!                    ***  ROUTINE asm_bgc_init_bkg  ***
217      !!
218      !! ** Purpose :   initialise bgc background
219      !!
220      !! ** Method  :   allocate arrays and read background from file
221      !!
222      !! ** Action  :   allocate arrays and read background from file
223      !!
224      !! References :   asm_inc_init
225      !!---------------------------------------------------------------------------
226      !!
227      INTEGER :: inum   ! i/o unit of background file
228      INTEGER :: jt     ! loop counter
229      !!
230      !!---------------------------------------------------------------------------
231
232#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
233      IF ( ln_slchltotinc ) THEN
234
235         ALLOCATE( pgrow_avg_bkg(jpi,jpj)        )
236         ALLOCATE( ploss_avg_bkg(jpi,jpj)        )
237         ALLOCATE( phyt_avg_bkg(jpi,jpj)         )
238         ALLOCATE( mld_max_bkg(jpi,jpj)          )
239         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
240         pgrow_avg_bkg(:,:)  = 0.0
241         ploss_avg_bkg(:,:)  = 0.0
242         phyt_avg_bkg(:,:)   = 0.0
243         mld_max_bkg(:,:)    = 0.0
244         tracer_bkg(:,:,:,:) = 0.0
245
246#if defined key_hadocc
247         ALLOCATE( chl_bkg(jpi,jpj)    )
248         ALLOCATE( cchl_p_bkg(jpi,jpj) )
249         chl_bkg(:,:)    = 0.0
250         cchl_p_bkg(:,:) = 0.0
251#endif
252         
253         !--------------------------------------------------------------------
254         ! Read background variables for slchltot assimilation
255         ! Some only required if performing balancing
256         !--------------------------------------------------------------------
257
258         CALL iom_open( c_asmbkg, inum )
259
260#if defined key_hadocc
261         CALL iom_get( inum, jpdom_autoglo, 'hadocc_chl',  chl_bkg    )
262         CALL iom_get( inum, jpdom_autoglo, 'hadocc_cchl', cchl_p_bkg )
263         chl_bkg(:,:)    = chl_bkg(:,:)    * tmask(:,:,1)
264         cchl_p_bkg(:,:) = cchl_p_bkg(:,:) * tmask(:,:,1)
265#elif defined key_medusa
266         CALL iom_get( inum, jpdom_autoglo, 'medusa_chn', tracer_bkg(:,:,:,jpchn) )
267         CALL iom_get( inum, jpdom_autoglo, 'medusa_chd', tracer_bkg(:,:,:,jpchd) )
268#endif
269         
270         IF ( ln_slchltotbal ) THEN
271
272            CALL iom_get( inum, jpdom_autoglo, 'pgrow_avg', pgrow_avg_bkg )
273            CALL iom_get( inum, jpdom_autoglo, 'ploss_avg', ploss_avg_bkg )
274            CALL iom_get( inum, jpdom_autoglo, 'phyt_avg',  phyt_avg_bkg  )
275            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
276            pgrow_avg_bkg(:,:) = pgrow_avg_bkg(:,:) * tmask(:,:,1)
277            ploss_avg_bkg(:,:) = ploss_avg_bkg(:,:) * tmask(:,:,1)
278            phyt_avg_bkg(:,:)  = phyt_avg_bkg(:,:)  * tmask(:,:,1)
279            mld_max_bkg(:,:)   = mld_max_bkg(:,:)   * tmask(:,:,1)
280
281#if defined key_hadocc
282            CALL iom_get( inum, jpdom_autoglo, 'hadocc_nut', tracer_bkg(:,:,:,jp_had_nut) )
283            CALL iom_get( inum, jpdom_autoglo, 'hadocc_phy', tracer_bkg(:,:,:,jp_had_phy) )
284            CALL iom_get( inum, jpdom_autoglo, 'hadocc_zoo', tracer_bkg(:,:,:,jp_had_zoo) )
285            CALL iom_get( inum, jpdom_autoglo, 'hadocc_pdn', tracer_bkg(:,:,:,jp_had_pdn) )
286            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
287            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
288#elif defined key_medusa
289            CALL iom_get( inum, jpdom_autoglo, 'medusa_phn', tracer_bkg(:,:,:,jpphn) )
290            CALL iom_get( inum, jpdom_autoglo, 'medusa_phd', tracer_bkg(:,:,:,jpphd) )
291            CALL iom_get( inum, jpdom_autoglo, 'medusa_pds', tracer_bkg(:,:,:,jppds) )
292            CALL iom_get( inum, jpdom_autoglo, 'medusa_zmi', tracer_bkg(:,:,:,jpzmi) )
293            CALL iom_get( inum, jpdom_autoglo, 'medusa_zme', tracer_bkg(:,:,:,jpzme) )
294            CALL iom_get( inum, jpdom_autoglo, 'medusa_din', tracer_bkg(:,:,:,jpdin) )
295            CALL iom_get( inum, jpdom_autoglo, 'medusa_sil', tracer_bkg(:,:,:,jpsil) )
296            CALL iom_get( inum, jpdom_autoglo, 'medusa_fer', tracer_bkg(:,:,:,jpfer) )
297            CALL iom_get( inum, jpdom_autoglo, 'medusa_det', tracer_bkg(:,:,:,jpdet) )
298            CALL iom_get( inum, jpdom_autoglo, 'medusa_dtc', tracer_bkg(:,:,:,jpdtc) )
299            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
300            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
301            CALL iom_get( inum, jpdom_autoglo, 'medusa_oxy', tracer_bkg(:,:,:,jpoxy) )
302#endif
303         ELSE IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
304#if defined key_hadocc
305            CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
306            CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
307#elif defined key_medusa
308            CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
309            CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
310#endif
311            CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
312            mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
313         ENDIF
314
315         CALL iom_close( inum )
316         
317         DO jt = 1, jptra
318            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
319         END DO
320     
321      ELSE IF ( ln_spco2inc .OR. ln_sfco2inc ) THEN
322
323         ALLOCATE( tracer_bkg(jpi,jpj,jpk,jptra) )
324         ALLOCATE( mld_max_bkg(jpi,jpj)          )
325         tracer_bkg(:,:,:,:) = 0.0
326         mld_max_bkg(:,:)    = 0.0
327
328         CALL iom_open( c_asmbkg, inum )
329         
330#if defined key_hadocc
331         CALL iom_get( inum, jpdom_autoglo, 'hadocc_dic', tracer_bkg(:,:,:,jp_had_dic) )
332         CALL iom_get( inum, jpdom_autoglo, 'hadocc_alk', tracer_bkg(:,:,:,jp_had_alk) )
333#elif defined key_medusa
334         CALL iom_get( inum, jpdom_autoglo, 'medusa_dic', tracer_bkg(:,:,:,jpdic) )
335         CALL iom_get( inum, jpdom_autoglo, 'medusa_alk', tracer_bkg(:,:,:,jpalk) )
336#endif
337         CALL iom_get( inum, jpdom_autoglo, 'mld_max',   mld_max_bkg   )
338
339         CALL iom_close( inum )
340         
341         DO jt = 1, jptra
342            tracer_bkg(:,:,:,jt) = tracer_bkg(:,:,:,jt) * tmask(:,:,:)
343         END DO
344         mld_max_bkg(:,:) = mld_max_bkg(:,:) * tmask(:,:,1)
345     
346      ENDIF
347#endif
348
349   END SUBROUTINE asm_bgc_init_bkg
350
351   !!======================================================================
352   !!======================================================================
353   !!======================================================================
354
355   SUBROUTINE asm_bgc_bal_wri( kt )
356      !!-----------------------------------------------------------------------
357      !!
358      !!                  ***  ROUTINE asm_bgc_bal_wri ***
359      !!
360      !! ** Purpose : Write to file the balancing increments
361      !!              calculated for biogeochemistry
362      !!
363      !! ** Method  : Write to file the balancing increments
364      !!              calculated for biogeochemistry
365      !!
366      !! ** Action  :
367      !!                   
368      !! References :
369      !!
370      !! History    :
371      !!        ! 2014-08 (D. Ford)  Initial version, based on asm_bkg_wri
372      !!-----------------------------------------------------------------------
373      !! * Arguments
374      INTEGER, INTENT( IN ) :: kt        ! Current time-step
375      !! * Local declarations
376      CHARACTER(LEN=50) :: cl_asmbal     ! Filename (with extension)
377      LOGICAL           :: llok          ! Check if file exists
378      INTEGER           :: inum          ! File unit number
379      REAL(wp)          :: zdate         ! Date
380      !!-----------------------------------------------------------------------
381     
382      ! Set things up
383      zdate = REAL( ndastp )
384      WRITE(cl_asmbal, FMT='(A,".nc")' ) TRIM( c_asmbal )
385
386      ! Check if file exists
387      INQUIRE( FILE = TRIM( cl_asmbal ), EXIST = llok )
388      IF( .NOT. llok ) THEN
389         IF(lwp) WRITE(numout,*) ' Setting up assimilation balancing increments file '// &
390            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
391
392         ! Define the output file       
393         CALL iom_open( c_asmbal, inum, ldwrt = .TRUE., kiolib = jprstlib)
394
395         ! Write the information
396         CALL iom_rstput( kt, kt, inum, 'rdastp' , zdate   )
397
398         IF ( ln_slchltotinc ) THEN
399#if defined key_medusa
400            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chn', slchltot_balinc(:,:,:,jpchn) )
401            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_chd', slchltot_balinc(:,:,:,jpchd) )
402            IF ( ln_slchltotbal ) THEN
403               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phn', slchltot_balinc(:,:,:,jpphn) )
404               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phd', slchltot_balinc(:,:,:,jpphd) )
405               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_pds', slchltot_balinc(:,:,:,jppds) )
406               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zmi', slchltot_balinc(:,:,:,jpzmi) )
407               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zme', slchltot_balinc(:,:,:,jpzme) )
408               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_din', slchltot_balinc(:,:,:,jpdin) )
409               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_sil', slchltot_balinc(:,:,:,jpsil) )
410               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_fer', slchltot_balinc(:,:,:,jpfer) )
411               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jpdet) )
412               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dtc', slchltot_balinc(:,:,:,jpdtc) )
413               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jpdic) )
414               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jpalk) )
415               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_oxy', slchltot_balinc(:,:,:,jpoxy) )
416            ENDIF
417#elif defined key_hadocc
418            CALL iom_rstput( kt, kt, inum, 'logchl_balinc_phy', slchltot_balinc(:,:,:,jp_had_phy) )
419            IF ( ln_slchltotbal ) THEN
420               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_nut', slchltot_balinc(:,:,:,jp_had_nut) )
421               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_zoo', slchltot_balinc(:,:,:,jp_had_zoo) )
422               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_det', slchltot_balinc(:,:,:,jp_had_pdn) )
423               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_dic', slchltot_balinc(:,:,:,jp_had_dic) )
424               CALL iom_rstput( kt, kt, inum, 'logchl_balinc_alk', slchltot_balinc(:,:,:,jp_had_alk) )
425            ENDIF
426#endif
427         ENDIF
428
429         IF ( ln_spco2inc ) THEN
430#if defined key_medusa
431            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jpdic) )
432            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jpalk) )
433#elif defined key_hadocc
434            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) )
435            CALL iom_rstput( kt, kt, inum, 'pco2_balinc_alk', spco2_balinc(:,:,:,jp_had_alk) )
436#endif
437         ELSE IF ( ln_sfco2inc ) THEN
438#if defined key_medusa
439            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jpdic) )
440            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_balinc(:,:,:,jpalk) )
441#elif defined key_hadocc
442            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_dic', spco2_balinc(:,:,:,jp_had_dic) )
443            CALL iom_rstput( kt, kt, inum, 'fco2_balinc_alk', spco2_balinc(:,:,:,jp_had_alk) )
444#endif
445         ENDIF
446
447         CALL iom_close( inum )
448      ELSE
449         CALL ctl_warn( TRIM( cl_asmbal ) // ' already exists ', &
450            &           ' Therefore not writing out balancing increments at this timestep', &
451            &           ' Something has probably gone wrong somewhere' )
452         IF(lwp) WRITE(numout,*) ' Failed to set up assimilation balancing increments file '// &
453            &                    TRIM( c_asmbal ) // ' at timestep = ', kt
454      ENDIF
455                                 
456   END SUBROUTINE asm_bgc_bal_wri
457
458   !!======================================================================
459   !!======================================================================
460   !!======================================================================
461
462   SUBROUTINE asm_bgc_bkg_wri( kt, knum )
463      !!---------------------------------------------------------------------------
464      !!                    ***  ROUTINE asm_bgc_bkg_wri  ***
465      !!
466      !! ** Purpose :   write out bgc background
467      !!
468      !! ** Method  :   write out bgc background
469      !!
470      !! ** Action  :   write out bgc background
471      !!
472      !! References :   asm_bkg_wri
473      !!---------------------------------------------------------------------------
474      !!
475      INTEGER, INTENT(in   ) :: kt     ! Current time-step
476      INTEGER, INTENT(in   ) :: knum   ! i/o unit of increments file
477      !!
478      !!---------------------------------------------------------------------------
479
480#if defined key_hadocc
481      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg             )
482      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg             )
483      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg              )
484      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max               )
485      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_nut'  , trn(:,:,:,jp_had_nut) )
486      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_phy'  , trn(:,:,:,jp_had_phy) )
487      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_zoo'  , trn(:,:,:,jp_had_zoo) )
488      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_pdn'  , trn(:,:,:,jp_had_pdn) )
489      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_dic'  , trn(:,:,:,jp_had_dic) )
490      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_alk'  , trn(:,:,:,jp_had_alk) )
491      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_chl'  , HADOCC_CHL(:,:,1)     )
492      CALL iom_rstput( kt, nitbkg_r, knum, 'hadocc_cchl' , cchl_p(:,:,1)         )
493#elif defined key_medusa && defined key_foam_medusa
494      CALL iom_rstput( kt, nitbkg_r, knum, 'pgrow_avg'   , pgrow_avg        )
495      CALL iom_rstput( kt, nitbkg_r, knum, 'ploss_avg'   , ploss_avg        )
496      CALL iom_rstput( kt, nitbkg_r, knum, 'phyt_avg'    , phyt_avg         )
497      CALL iom_rstput( kt, nitbkg_r, knum, 'mld_max'     , mld_max          )
498      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chn'  , trn(:,:,:,jpchn) )
499      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_chd'  , trn(:,:,:,jpchd) )
500      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phn'  , trn(:,:,:,jpphn) )
501      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_phd'  , trn(:,:,:,jpphd) )
502      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_pds'  , trn(:,:,:,jppds) )
503      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zmi'  , trn(:,:,:,jpzmi) )
504      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_zme'  , trn(:,:,:,jpzme) )
505      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_din'  , trn(:,:,:,jpdin) )
506      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_sil'  , trn(:,:,:,jpsil) )
507      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_fer'  , trn(:,:,:,jpfer) )
508      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_det'  , trn(:,:,:,jpdet) )
509      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dtc'  , trn(:,:,:,jpdtc) )
510      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_dic'  , trn(:,:,:,jpdic) )
511      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_alk'  , trn(:,:,:,jpalk) )
512      CALL iom_rstput( kt, nitbkg_r, knum, 'medusa_oxy'  , trn(:,:,:,jpoxy) )
513#endif
514
515   END SUBROUTINE asm_bgc_bkg_wri
516
517   !!======================================================================
518   !!======================================================================
519   !!======================================================================
520
521   SUBROUTINE slchltot_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau )
522      !!----------------------------------------------------------------------
523      !!                    ***  ROUTINE logchl_asm_inc  ***
524      !!         
525      !! ** Purpose : Apply the chlorophyll assimilation increments.
526      !!
527      !! ** Method  : Calculate increments to state variables using nitrogen
528      !!              balancing.
529      !!              Direct initialization or Incremental Analysis Updating.
530      !!
531      !! ** Action  :
532      !!----------------------------------------------------------------------
533      INTEGER,  INTENT(IN) :: kt        ! Current time step
534      LOGICAL,  INTENT(IN) :: ll_asmdin ! Flag for direct initialisation
535      LOGICAL,  INTENT(IN) :: ll_asmiau ! Flag for incremental analysis update
536      INTEGER,  INTENT(IN) :: kcycper   ! Dimension of pwgtiau
537      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
538      !
539      INTEGER  :: jk              ! Loop counter
540      INTEGER  :: it              ! Index
541      REAL(wp) :: zincwgt         ! IAU weight for current time step
542      REAL(wp) :: zincper         ! IAU interval in seconds
543      !!----------------------------------------------------------------------
544
545      IF ( kt <= nit000 ) THEN
546
547         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt
548
549#if defined key_medusa && defined key_foam_medusa
550         CALL asm_logchl_bal_medusa( slchltot_bkginc, zincper, mld_choice_bgc, &
551            &                        rn_maxchlinc, ln_slchltotbal, ll_asmdin,  &
552            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
553            &                        phyt_avg_bkg, mld_max_bkg,              &
554            &                        tracer_bkg, slchltot_balinc )
555#elif defined key_hadocc
556         CALL asm_logchl_bal_hadocc( slchltot_bkginc, zincper, mld_choice_bgc, &
557            &                        rn_maxchlinc, ln_slchltotbal, ll_asmdin,  &
558            &                        pgrow_avg_bkg, ploss_avg_bkg,           &
559            &                        phyt_avg_bkg, mld_max_bkg,              &
560            &                        chl_bkg, cchl_p_bkg,                    &
561            &                        tracer_bkg, slchltot_balinc )
562#else
563         CALL ctl_stop( 'Attempting to assimilate slchltot, ', &
564            &           'but not defined a biogeochemical model' )
565#endif
566
567      ENDIF
568
569      IF ( ll_asmiau ) THEN
570
571         !--------------------------------------------------------------------
572         ! Incremental Analysis Updating
573         !--------------------------------------------------------------------
574
575         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
576
577            it = kt - nit000 + 1
578            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
579            ! note this is not a tendency so should not be divided by rdt
580
581            IF(lwp) THEN
582               WRITE(numout,*) 
583               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', &
584                  &  kt,' with IAU weight = ', pwgtiau(it)
585               WRITE(numout,*) '~~~~~~~~~~~~'
586            ENDIF
587
588            ! Update the biogeochemical variables
589            ! Add directly to trn and trb, rather than to tra, because tra gets
590            ! reset to zero at the start of trc_stp, called after this routine
591#if defined key_medusa && defined key_foam_medusa
592            DO jk = 1, jpkm1
593               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
594                  &                          slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
595               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
596                  &                          slchltot_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
597            END DO
598#elif defined key_hadocc
599            DO jk = 1, jpkm1
600               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
601                  &                          slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
602               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
603                  &                          slchltot_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
604            END DO
605#endif
606
607            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
608            ! which is called at end of model run
609         ENDIF
610
611      ELSEIF ( ll_asmdin ) THEN 
612
613         !--------------------------------------------------------------------
614         ! Direct Initialization
615         !--------------------------------------------------------------------
616         
617         IF ( kt == nitdin_r ) THEN
618
619            neuler = 0                    ! Force Euler forward step
620
621#if defined key_medusa && defined key_foam_medusa
622            ! Initialize the now fields with the background + increment
623            ! Background currently is what the model is initialised with
624            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', &
625               &           ' Background state is taken from model rather than background file' )
626            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
627               &                         slchltot_balinc(:,:,:,jp_msa0:jp_msa1)
628            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
629#elif defined key_hadocc
630            ! Initialize the now fields with the background + increment
631            ! Background currently is what the model is initialised with
632            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', &
633               &           ' Background state is taken from model rather than background file' )
634            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
635               &                         slchltot_balinc(:,:,:,jp_had0:jp_had1)
636            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
637#endif
638 
639            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
640            ! which is called at end of model run
641         ENDIF
642         !
643      ENDIF
644      !
645   END SUBROUTINE slchltot_asm_inc
646
647   !!======================================================================
648   !!======================================================================
649   !!======================================================================
650
651   SUBROUTINE spco2_asm_inc( kt, ll_asmdin, ll_asmiau, kcycper, pwgtiau, &
652      &                      ll_trainc, pt_bkginc, ps_bkginc )
653      !!----------------------------------------------------------------------
654      !!                    ***  ROUTINE pco2_asm_inc  ***
655      !!         
656      !! ** Purpose : Apply the pco2/fco2 assimilation increments.
657      !!
658      !! ** Method  : Calculate increments to state variables using carbon
659      !!              balancing.
660      !!              Direct initialization or Incremental Analysis Updating.
661      !!
662      !! ** Action  :
663      !!----------------------------------------------------------------------
664      INTEGER, INTENT(IN)                   :: kt           ! Current time step
665      LOGICAL, INTENT(IN)                   :: ll_asmdin    ! Flag for direct initialisation
666      LOGICAL, INTENT(IN)                   :: ll_asmiau    ! Flag for incremental analysis update
667      INTEGER, INTENT(IN)                   :: kcycper      ! Dimension of pwgtiau
668      REAL(wp), DIMENSION(kcycper), INTENT(IN) :: pwgtiau   ! IAU weights
669      LOGICAL, INTENT(IN)                   :: ll_trainc    ! Flag for T/S increments
670      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: pt_bkginc ! T increments
671      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ps_bkginc ! S increments
672      !
673      INTEGER                               :: ji, jj, jk   ! Loop counters
674      INTEGER                               :: it           ! Index
675      INTEGER                               :: jkmax        ! Index of mixed layer depth
676      REAL(wp)                              :: zincwgt      ! IAU weight for current time step
677      REAL(wp), DIMENSION(jpi,jpj)          :: zmld         ! Mixed layer depth for balancing
678      REAL(wp), DIMENSION(jpi,jpj)          :: pco2_bkginc_temp ! For fCO2 conversion
679      REAL(wp), DIMENSION(jpi,jpj)          :: dic_bkg_temp ! DIC background
680      REAL(wp), DIMENSION(jpi,jpj)          :: alk_bkg_temp ! Alkalinity background
681      REAL(wp), DIMENSION(jpi,jpj)          :: tem_bkg_temp ! Temperature background
682      REAL(wp), DIMENSION(jpi,jpj)          :: sal_bkg_temp ! Salinity background
683      REAL(wp), DIMENSION(1)                :: patm         ! Atmospheric pressure
684      REAL(wp), DIMENSION(1)                :: phyd         ! Hydrostatic pressure
685
686      ! Coefficients for fCO2 to pCO2 conversion
687      REAL(wp)                              :: zcoef_fco2_1 = -1636.75
688      REAL(wp)                              :: zcoef_fco2_2 = 12.0408
689      REAL(wp)                              :: zcoef_fco2_3 = 0.0327957
690      REAL(wp)                              :: zcoef_fco2_4 = 0.0000316528
691      REAL(wp)                              :: zcoef_fco2_5 = 2.0
692      REAL(wp)                              :: zcoef_fco2_6 = 57.7
693      REAL(wp)                              :: zcoef_fco2_7 = 0.118
694      REAL(wp)                              :: zcoef_fco2_8 = 82.0578
695      !!----------------------------------------------------------------------
696
697      IF ( kt <= nit000 ) THEN
698
699         !--------------------------------------------------------------------
700         ! DIC and alkalinity balancing
701         !--------------------------------------------------------------------
702
703         IF ( ln_sfco2inc ) THEN
704#if defined key_medusa && defined key_roam
705            ! If assimilating fCO2, then convert to pCO2 using MEDUSA MOCSY subroutine
706            patm(1) = 1.0
707            phyd(1) = 0.0
708            DO jj = 1, jpj
709               DO ji = 1, jpi
710                  CALL f2pCO2( spco2_bkginc(ji,jj), tsn(ji,jj,1,1), patm, phyd, 1, pco2_bkginc_temp(ji,jj) )
711               END DO
712            END DO
713#else
714            ! If assimilating fCO2, then convert to pCO2 using temperature
715            ! See flux_gas.F90 within HadOCC for details of calculation
716            pco2_bkginc_temp(:,:) = spco2_bkginc(:,:) /                                                             &
717               &                    EXP((zcoef_fco2_1                                                            + &
718               &                         zcoef_fco2_2 * (tsn(:,:,1,1)+rt0)                                       - &
719               &                         zcoef_fco2_3 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)                    + &
720               &                         zcoef_fco2_4 * (tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0)*(tsn(:,:,1,1)+rt0) + &
721               &                         zcoef_fco2_5 * (zcoef_fco2_6 - zcoef_fco2_7 * (tsn(:,:,1,1)+rt0)))      / &
722               &                        (zcoef_fco2_8 * (tsn(:,:,1,1)+rt0)))
723#endif
724         ELSE
725            pco2_bkginc_temp(:,:) = spco2_bkginc(:,:)
726         ENDIF
727
728         ! Account for physics assimilation if required
729         IF ( ll_trainc ) THEN
730            tem_bkg_temp(:,:) = tsn(:,:,1,1) + pt_bkginc(:,:,1)
731            sal_bkg_temp(:,:) = tsn(:,:,1,2) + ps_bkginc(:,:,1)
732         ELSE
733            tem_bkg_temp(:,:) = tsn(:,:,1,1)
734            sal_bkg_temp(:,:) = tsn(:,:,1,2)
735         ENDIF
736
737#if defined key_medusa
738         ! Account for logchl balancing if required
739         IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN
740            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic) + slchltot_balinc(:,:,1,jpdic)
741            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk) + slchltot_balinc(:,:,1,jpalk)
742         ELSE
743            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jpdic)
744            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jpalk)
745         ENDIF
746
747         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
748            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
749            &               spco2_balinc(:,:,1,jpdic), spco2_balinc(:,:,1,jpalk) )
750
751#elif defined key_hadocc
752         ! Account for slchltot balancing if required
753         IF ( ln_slchltotinc .AND. ln_slchltotbal ) THEN
754            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic) + slchltot_balinc(:,:,1,jp_had_dic)
755            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk) + slchltot_balinc(:,:,1,jp_had_alk)
756         ELSE
757            dic_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_dic)
758            alk_bkg_temp(:,:) = tracer_bkg(:,:,1,jp_had_alk)
759         ENDIF
760
761         CALL asm_pco2_bal( pco2_bkginc_temp(:,:), dic_bkg_temp(:,:), alk_bkg_temp(:,:), &
762            &               tem_bkg_temp(:,:), sal_bkg_temp(:,:),                        &
763            &               spco2_balinc(:,:,1,jp_had_dic), spco2_balinc(:,:,1,jp_had_alk) )
764
765#else
766         CALL ctl_stop( 'Attempting to assimilate pCO2/fCO2, ', &
767            &           'but not defined a biogeochemical model' )
768#endif
769
770         ! Select mixed layer
771         IF ( ll_asmdin ) THEN
772#if defined key_hadocc || (defined key_medusa && defined key_foam_medusa)
773            CALL ctl_warn( ' Doing direct initialisation with pCO2 assimilation', &
774               &           ' Mixed layer depth taken to be background maximum mld_max_bkg' )
775            zmld(:,:) = mld_max_bkg(:,:)
776#else
777            CALL ctl_stop( 'No biogeochemical model defined for pCO2 assimilation' )
778#endif
779         ELSE
780            SELECT CASE( mld_choice_bgc )
781            CASE ( 1 )                   ! Turbocline/mixing depth [W points]
782               zmld(:,:) = hmld(:,:)
783            CASE ( 2 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [W points]
784               zmld(:,:) = hmlp(:,:)
785            CASE ( 3 )                   ! Kara MLD [Interpolated]
786#if defined key_karaml
787               IF ( ln_kara ) THEN
788                  zmld(:,:) = hmld_kara(:,:)
789               ELSE
790                  CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
791                     &           ' but ln_kara=.false.' )
792               ENDIF
793#else
794               CALL ctl_stop( ' Kara mixed layer requested for pCO2 assimilation,', &
795                  &           ' but is not defined' )
796#endif
797            CASE ( 4 )                   ! Temperature criterion (0.2 K change from surface) [T points]
798               !zmld(:,:) = hmld_tref(:,:)
799               CALL ctl_stop( ' hmld_tref mixed layer requested for pCO2 assimilation,', &
800                  &           ' but is not available in this version' )
801            CASE ( 5 )                   ! Density criterion (0.01 kg/m^3 change from 10m) [T points]
802               zmld(:,:) = hmlpt(:,:)
803            END SELECT
804         ENDIF
805
806         ! Propagate through mixed layer
807         DO jj = 1, jpj
808            DO ji = 1, jpi
809               !
810               jkmax = jpk-1
811               DO jk = jpk-1, 1, -1
812                  IF ( ( zmld(ji,jj) >  gdepw_n(ji,jj,jk)   ) .AND. &
813                     & ( zmld(ji,jj) <= gdepw_n(ji,jj,jk+1) ) ) THEN
814                     zmld(ji,jj) = gdepw_n(ji,jj,jk+1)
815                     jkmax = jk
816                  ENDIF
817               END DO
818               !
819#if defined key_top
820               DO jk = 2, jkmax
821                  spco2_balinc(ji,jj,jk,:) = spco2_balinc(ji,jj,1,:)
822               END DO
823#endif
824               !
825            END DO
826         END DO
827
828      ENDIF
829
830      IF ( ll_asmiau ) THEN
831
832         !--------------------------------------------------------------------
833         ! Incremental Analysis Updating
834         !--------------------------------------------------------------------
835
836         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
837
838            it = kt - nit000 + 1
839            zincwgt = pwgtiau(it)      ! IAU weight for the current time step
840            ! note this is not a tendency so should not be divided by rdt
841
842            IF(lwp) THEN
843               WRITE(numout,*) 
844               IF ( ln_spco2inc ) THEN
845                  WRITE(numout,*) 'pco2_asm_inc : pco2 IAU at time step = ', &
846                     &  kt,' with IAU weight = ', pwgtiau(it)
847               ELSE IF ( ln_sfco2inc ) THEN
848                  WRITE(numout,*) 'fco2_asm_inc : fco2 IAU at time step = ', &
849                     &  kt,' with IAU weight = ', pwgtiau(it)
850               ENDIF
851               WRITE(numout,*) '~~~~~~~~~~~~'
852            ENDIF
853
854            ! Update the biogeochemical variables
855            ! Add directly to trn and trb, rather than to tra, because tra gets
856            ! reset to zero at the start of trc_stp, called after this routine
857#if defined key_medusa && defined key_foam_medusa
858            DO jk = 1, jpkm1
859               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + &
860                  &                          spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
861               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + &
862                  &                          spco2_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt
863            END DO
864#elif defined key_hadocc
865            DO jk = 1, jpkm1
866               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + &
867                  &                          spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
868               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + &
869                  &                          spco2_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt
870            END DO
871#endif
872
873            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
874            ! which is called at end of model run
875
876         ENDIF
877
878      ELSEIF ( ll_asmdin ) THEN 
879
880         !--------------------------------------------------------------------
881         ! Direct Initialization
882         !--------------------------------------------------------------------
883         
884         IF ( kt == nitdin_r ) THEN
885
886            neuler = 0                    ! Force Euler forward step
887
888#if defined key_medusa && defined key_foam_medusa
889            ! Initialize the now fields with the background + increment
890            ! Background currently is what the model is initialised with
891            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with pCO2 assimilation', &
892               &           ' Background state is taken from model rather than background file' )
893            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + &
894               &                         spco2_balinc(:,:,:,jp_msa0:jp_msa1)
895            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1)
896#elif defined key_hadocc
897            ! Initialize the now fields with the background + increment
898            ! Background currently is what the model is initialised with
899            CALL ctl_warn( ' Doing direct initialisation of HadOCC with pCO2 assimilation', &
900               &           ' Background state is taken from model rather than background file' )
901            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + &
902               &                         spco2_balinc(:,:,:,jp_had0:jp_had1)
903            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1)
904#endif
905 
906            ! Do not deallocate arrays - needed by asm_bgc_bal_wri
907            ! which is called at end of model run
908         ENDIF
909         !
910      ENDIF
911      !
912   END SUBROUTINE spco2_asm_inc
913
914   !!======================================================================
915
916END MODULE asmbgc
Note: See TracBrowser for help on using the repository browser.