source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 10149

Last change on this file since 10149 was 10149, checked in by frrh, 23 months ago

Met Office GMED ticket 379: Merged David Ford's MEDUSA assimilation changes
using command:

svn merge -r 10054:10141 svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package_asm_3d_bgc_v3

File size: 55.1 KB
Line 
1MODULE asminc
2   !!======================================================================
3   !!                       ***  MODULE asminc  ***
4   !! Assimilation increment : Apply an increment generated by data
5   !!                          assimilation
6   !!======================================================================
7   !! History :       ! 2007-03  (M. Martin)  Met Office version
8   !!                 ! 2007-04  (A. Weaver)  calc_date original code
9   !!                 ! 2007-04  (A. Weaver)  Merge with OPAVAR/NEMOVAR
10   !!   NEMO     3.3  ! 2010-05  (D. Lea)  Update to work with NEMO v3.2
11   !!             -   ! 2010-05  (D. Lea)  add calc_month_len routine based on day_init
12   !!            3.4  ! 2012-10  (A. Weaver and K. Mogensen) Fix for direct initialization
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   'key_asminc'   : Switch on the assimilation increment interface
17   !!----------------------------------------------------------------------
18   !!   asm_inc_init   : Initialize the increment arrays and IAU weights
19   !!   calc_date      : Compute the calendar date YYYYMMDD on a given step
20   !!   tra_asm_inc    : Apply the tracer (T and S) increments
21   !!   dyn_asm_inc    : Apply the dynamic (u and v) increments
22   !!   ssh_asm_inc    : Apply the SSH increment
23   !!   seaice_asm_inc : Apply the seaice increment
24   !!----------------------------------------------------------------------
25   USE wrk_nemo         ! Memory Allocation
26   USE par_oce          ! Ocean space and time domain variables
27   USE dom_oce          ! Ocean space and time domain
28   USE domvvl           ! domain: variable volume level
29   USE oce              ! Dynamics and active tracers defined in memory
30   USE ldfdyn_oce       ! ocean dynamics: lateral physics
31   USE eosbn2           ! Equation of state - in situ and potential density
32   USE zpshde           ! Partial step : Horizontal Derivative
33   USE iom              ! Library to read input files
34   USE asmpar           ! Parameters for the assmilation interface
35   USE c1d              ! 1D initialization
36   USE in_out_manager   ! I/O manager
37   USE lib_mpp          ! MPP library
38#if defined key_lim2
39   USE ice_2            ! LIM2
40#endif
41#if defined key_cice && defined key_asminc
42   USE sbc_ice, ONLY : & ! CICE Ice model variables
43   & ndaice_da, nfresh_da, nfsalt_da
44#endif
45   USE sbc_oce          ! Surface boundary condition variables.
46   USE asmbgc           ! Biogeochemistry assimilation
47
48   IMPLICIT NONE
49   PRIVATE
50   
51   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
52   PUBLIC   calc_date      !: Compute the calendar date YYYYMMDD on a given step
53   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
54   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
55   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
56   PUBLIC   seaice_asm_inc !: Apply the seaice increment
57   PUBLIC   bgc_asm_inc    !: Apply the biogeochemistry increments
58
59#if defined key_asminc
60    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
61#else
62    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
63#endif
64   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields
65   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment
66   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization
67   LOGICAL, PUBLIC :: ln_trainc = .FALSE.      !: No tracer (T and S) assimilation increments
68   LOGICAL, PUBLIC :: ln_dyninc = .FALSE.      !: No dynamics (u and v) assimilation increments
69   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment
70   LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment
71   LOGICAL, PUBLIC :: lk_bgcinc = .FALSE.      !: No biogeochemistry increments
72   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check
73   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
74   INTEGER, PUBLIC :: nn_divdmp                !: Apply divergence damping filter nn_divdmp times
75
76   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
77   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
78   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
79   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
80   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
81#if defined key_asminc
82   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_iau           !: IAU-weighted sea surface height increment
83#endif
84   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
85   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
86   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
87   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
88   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
89   !
90   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
91   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
92   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
93
94   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
95   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
96
97   !! * Substitutions
98#  include "domzgr_substitute.h90"
99#  include "ldfdyn_substitute.h90"
100#  include "vectopt_loop_substitute.h90"
101   !!----------------------------------------------------------------------
102   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
103   !! $Id$
104   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
105   !!----------------------------------------------------------------------
106CONTAINS
107
108   SUBROUTINE asm_inc_init
109      !!----------------------------------------------------------------------
110      !!                    ***  ROUTINE asm_inc_init  ***
111      !!         
112      !! ** Purpose : Initialize the assimilation increment and IAU weights.
113      !!
114      !! ** Method  : Initialize the assimilation increment and IAU weights.
115      !!
116      !! ** Action  :
117      !!----------------------------------------------------------------------
118      INTEGER :: ji, jj, jk, jt  ! dummy loop indices
119      INTEGER :: imid, inum      ! local integers
120      INTEGER :: ios             ! Local integer output status for namelist read
121      INTEGER :: iiauper         ! Number of time steps in the IAU period
122      INTEGER :: icycper         ! Number of time steps in the cycle
123      INTEGER :: iitend_date     ! Date YYYYMMDD of final time step
124      INTEGER :: iitbkg_date     ! Date YYYYMMDD of background time step for Jb term
125      INTEGER :: iitdin_date     ! Date YYYYMMDD of background time step for DI
126      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step
127      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step
128      !
129      REAL(wp) :: znorm        ! Normalization factor for IAU weights
130      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights (should be equal to one)
131      REAL(wp) :: z_inc_dateb  ! Start date of interval on which increment is valid
132      REAL(wp) :: z_inc_datef  ! End date of interval on which increment is valid
133      REAL(wp) :: zdate_bkg    ! Date in background state file for DI
134      REAL(wp) :: zdate_inc    ! Time axis in increments file
135      !
136      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace
137      !!
138      NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri,                           &
139         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
140         &                 ln_phytobal, ln_slchltotinc, ln_slchldiainc,    &
141         &                 ln_slchlnoninc, ln_schltotinc, ln_slphytotinc,  &
142         &                 ln_slphydiainc, ln_slphynoninc, ln_spco2inc,    &
143         &                 ln_sfco2inc, ln_plchltotinc, ln_pchltotinc,     &
144         &                 ln_pno3inc, ln_psi4inc, ln_pdicinc, ln_palkinc, &
145         &                 ln_pphinc, ln_po2inc,                           &
146         &                 ln_asmdin, ln_asmiau,                           &
147         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
148         &                 ln_salfix, salfixmin, nn_divdmp,                &
149         &                 ln_seaiceinc, ln_temnofreeze,                   &
150         &                 mld_choice_bgc, rn_maxchlinc
151      !!----------------------------------------------------------------------
152
153      !-----------------------------------------------------------------------
154      ! Read Namelist nam_asminc : assimilation increment interface
155      !-----------------------------------------------------------------------
156      ln_seaiceinc = .FALSE.
157      ln_temnofreeze = .FALSE.
158
159      REWIND( numnam_ref )              ! Namelist nam_asminc in reference namelist : Assimilation increment
160      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
161901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in reference namelist', lwp )
162
163      REWIND( numnam_cfg )              ! Namelist nam_asminc in configuration namelist : Assimilation increment
164      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
165902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_asminc in configuration namelist', lwp )
166      IF(lwm) WRITE ( numond, nam_asminc )
167
168      ! Control print
169      IF(lwp) THEN
170         WRITE(numout,*)
171         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
172         WRITE(numout,*) '~~~~~~~~~~~~'
173         WRITE(numout,*) '   Namelist nam_asminc : set assimilation increment parameters'
174         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
175         WRITE(numout,*) '      Logical switch for writing out balancing increments      ln_balwri = ', ln_balwri
176         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
177         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
178         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
179         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
180         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
181         WRITE(numout,*) '      Logical switch for phytoplankton balancing             ln_phytobal = ', ln_phytobal
182         WRITE(numout,*) '      Logical switch for applying slchltot increments     ln_slchltotinc = ', ln_slchltotinc
183         WRITE(numout,*) '      Logical switch for applying slchldia increments     ln_slchldiainc = ', ln_slchldiainc
184         WRITE(numout,*) '      Logical switch for applying slchlnon increments     ln_slchlnoninc = ', ln_slchlnoninc
185         WRITE(numout,*) '      Logical switch for applying schltot increments       ln_schltotinc = ', ln_schltotinc
186         WRITE(numout,*) '      Logical switch for applying slphytot increments     ln_slphytotinc = ', ln_slphytotinc
187         WRITE(numout,*) '      Logical switch for applying slphydia increments     ln_slphydiainc = ', ln_slphydiainc
188         WRITE(numout,*) '      Logical switch for applying slphynon increments     ln_slphynoninc = ', ln_slphynoninc
189         WRITE(numout,*) '      Logical switch for applying spco2 increments           ln_spco2inc = ', ln_spco2inc
190         WRITE(numout,*) '      Logical switch for applying sfco2 increments           ln_sfco2inc = ', ln_sfco2inc
191         WRITE(numout,*) '      Logical switch for applying plchltot increments     ln_plchltotinc = ', ln_plchltotinc
192         WRITE(numout,*) '      Logical switch for applying pchltot increments       ln_pchltotinc = ', ln_pchltotinc
193         WRITE(numout,*) '      Logical switch for applying pno3 increments             ln_pno3inc = ', ln_pno3inc
194         WRITE(numout,*) '      Logical switch for applying psi4 increments             ln_psi4inc = ', ln_psi4inc
195         WRITE(numout,*) '      Logical switch for applying pdic increments             ln_pdicinc = ', ln_pdicinc
196         WRITE(numout,*) '      Logical switch for applying palk increments             ln_palkinc = ', ln_palkinc
197         WRITE(numout,*) '      Logical switch for applying pph increments               ln_pphinc = ', ln_pphinc
198         WRITE(numout,*) '      Logical switch for applying po2 increments               ln_po2inc = ', ln_po2inc
199         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
200         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
201         WRITE(numout,*) '      Timestep of background for DI in [0,nitend-nit000-1]     nitdin    = ', nitdin
202         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
203         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
204         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
205         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
206         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
207         WRITE(numout,*) '      Choice of MLD for BGC assimilation                  mld_choice_bgc = ', mld_choice_bgc
208         WRITE(numout,*) '      Maximum absolute chlorophyll increment (<=0 = off)    rn_maxchlinc = ', rn_maxchlinc
209      ENDIF
210
211      nitbkg_r    = nitbkg    + nit000 - 1  ! Background time referenced to nit000
212      nitdin_r    = nitdin    + nit000 - 1  ! Background time for DI referenced to nit000
213      nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000
214      nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000
215
216      iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length
217      icycper = nitend      - nit000      + 1  ! Cycle interval length
218
219      CALL calc_date( nit000, nitend     , ndate0, iitend_date    )     ! Date of final time step
220      CALL calc_date( nit000, nitbkg_r   , ndate0, iitbkg_date    )     ! Background time for Jb referenced to ndate0
221      CALL calc_date( nit000, nitdin_r   , ndate0, iitdin_date    )     ! Background time for DI referenced to ndate0
222      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )     ! IAU start time referenced to ndate0
223      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )     ! IAU end time referenced to ndate0
224      !
225      IF(lwp) THEN
226         WRITE(numout,*)
227         WRITE(numout,*) '   Time steps referenced to current cycle:'
228         WRITE(numout,*) '       iitrst      = ', nit000 - 1
229         WRITE(numout,*) '       nit000      = ', nit000
230         WRITE(numout,*) '       nitend      = ', nitend
231         WRITE(numout,*) '       nitbkg_r    = ', nitbkg_r
232         WRITE(numout,*) '       nitdin_r    = ', nitdin_r
233         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
234         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
235         WRITE(numout,*)
236         WRITE(numout,*) '   Dates referenced to current cycle:'
237         WRITE(numout,*) '       ndastp         = ', ndastp
238         WRITE(numout,*) '       ndate0         = ', ndate0
239         WRITE(numout,*) '       iitend_date    = ', iitend_date
240         WRITE(numout,*) '       iitbkg_date    = ', iitbkg_date
241         WRITE(numout,*) '       iitdin_date    = ', iitdin_date
242         WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date
243         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date
244      ENDIF
245     
246      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
247         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
248         & ln_slphynoninc .OR. ln_spco2inc    .OR. ln_sfco2inc    .OR. &
249         & ln_plchltotinc .OR. ln_pchltotinc  .OR. ln_pno3inc     .OR. &
250         & ln_psi4inc     .OR. ln_pdicinc     .OR. ln_palkinc     .OR. &
251         & ln_pphinc      .OR. ln_po2inc ) THEN
252         lk_bgcinc = .TRUE.
253      ENDIF
254
255      IF ( nacc /= 0 ) &
256         & CALL ctl_stop( ' nacc /= 0 and key_asminc :',  &
257         &                ' Assimilation increments have only been implemented', &
258         &                ' for synchronous time stepping' )
259
260      IF ( ( ln_asmdin ).AND.( ln_asmiau ) )   &
261         & CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
262         &                ' Choose Direct Initialization OR Incremental Analysis Updating')
263
264      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
265         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. &
266         &        ( lk_bgcinc ) )) &
267         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
268         &                ' ln_(bgc-variable)inc is set to .true.', &
269         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
270         &                ' Inconsistent options')
271
272      IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
273         & CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :',  &
274         &                ' Type IAU weighting function is invalid')
275
276      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
277         & .AND.( .NOT. lk_bgcinc ) )  &
278         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', &
279         &                ' ln_(bgc-variable)inc are set to .false. :', &
280         &                ' The assimilation increments are not applied')
281
282      IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
283         & CALL ctl_stop( ' nitiaustr = nitiaufin :',  &
284         &                ' IAU interval is of zero length')
285
286      IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
287         & CALL ctl_stop( ' nitiaustr or nitiaufin :',  &
288         &                ' IAU starting or final time step is outside the cycle interval', &
289         &                 ' Valid range nit000 to nitend')
290
291      IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
292         & CALL ctl_stop( ' nitbkg :',  &
293         &                ' Background time step is outside the cycle interval')
294
295      IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
296         & CALL ctl_stop( ' nitdin :',  &
297         &                ' Background time step for Direct Initialization is outside', &
298         &                ' the cycle interval')
299     
300      IF ( lk_bgcinc ) CALL asm_bgc_check_options
301
302      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
303
304      !--------------------------------------------------------------------
305      ! Initialize the Incremental Analysis Updating weighting function
306      !--------------------------------------------------------------------
307
308      IF ( ln_asmiau ) THEN
309
310         ALLOCATE( wgtiau( icycper ) )
311
312         wgtiau(:) = 0.0
313
314         IF ( niaufn == 0 ) THEN
315
316            !---------------------------------------------------------
317            ! Constant IAU forcing
318            !---------------------------------------------------------
319
320            DO jt = 1, iiauper
321               wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
322            END DO
323
324         ELSEIF ( niaufn == 1 ) THEN
325
326            !---------------------------------------------------------
327            ! Linear hat-like, centred in middle of IAU interval
328            !---------------------------------------------------------
329
330            ! Compute the normalization factor
331            znorm = 0.0
332            IF ( MOD( iiauper, 2 ) == 0 ) THEN  ! Even number of time steps in IAU interval
333               imid = iiauper / 2 
334               DO jt = 1, imid
335                  znorm = znorm + REAL( jt )
336               END DO
337               znorm = 2.0 * znorm
338            ELSE                               ! Odd number of time steps in IAU interval
339               imid = ( iiauper + 1 ) / 2       
340               DO jt = 1, imid - 1
341                  znorm = znorm + REAL( jt )
342               END DO
343               znorm = 2.0 * znorm + REAL( imid )
344            ENDIF
345            znorm = 1.0 / znorm
346
347            DO jt = 1, imid - 1
348               wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
349            END DO
350            DO jt = imid, iiauper
351               wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
352            END DO
353
354         ENDIF
355
356         ! Test that the integral of the weights over the weighting interval equals 1
357          IF(lwp) THEN
358             WRITE(numout,*)
359             WRITE(numout,*) 'asm_inc_init : IAU weights'
360             WRITE(numout,*) '~~~~~~~~~~~~'
361             WRITE(numout,*) '             time step         IAU  weight'
362             WRITE(numout,*) '             =========     ====================='
363             ztotwgt = 0.0
364             DO jt = 1, icycper
365                ztotwgt = ztotwgt + wgtiau(jt)
366                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
367             END DO   
368             WRITE(numout,*) '         ==================================='
369             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
370             WRITE(numout,*) '         ==================================='
371          ENDIF
372         
373      ENDIF
374
375      !--------------------------------------------------------------------
376      ! Allocate and initialize the increment arrays
377      !--------------------------------------------------------------------
378
379      ALLOCATE( t_bkginc(jpi,jpj,jpk) )
380      ALLOCATE( s_bkginc(jpi,jpj,jpk) )
381      ALLOCATE( u_bkginc(jpi,jpj,jpk) )
382      ALLOCATE( v_bkginc(jpi,jpj,jpk) )
383      ALLOCATE( ssh_bkginc(jpi,jpj)   )
384      ALLOCATE( seaice_bkginc(jpi,jpj))
385#if defined key_asminc
386      ALLOCATE( ssh_iau(jpi,jpj)      )
387#endif
388      t_bkginc(:,:,:) = 0.0
389      s_bkginc(:,:,:) = 0.0
390      u_bkginc(:,:,:) = 0.0
391      v_bkginc(:,:,:) = 0.0
392      ssh_bkginc(:,:) = 0.0
393      seaice_bkginc(:,:) = 0.0
394#if defined key_asminc
395      ssh_iau(:,:)    = 0.0
396#endif
397      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) &
398         &  .OR.( lk_bgcinc ) ) THEN
399
400         !--------------------------------------------------------------------
401         ! Read the increments from file
402         !--------------------------------------------------------------------
403
404         CALL iom_open( c_asminc, inum )
405
406         CALL iom_get( inum, 'time', zdate_inc ) 
407
408         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
409         CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
410         z_inc_dateb = zdate_inc
411         z_inc_datef = zdate_inc
412
413         IF(lwp) THEN
414            WRITE(numout,*) 
415            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
416               &            ' between dates ', NINT( z_inc_dateb ),' and ',  &
417               &            NINT( z_inc_datef )
418            WRITE(numout,*) '~~~~~~~~~~~~'
419         ENDIF
420
421         IF (     ( NINT( z_inc_dateb ) < ndastp      ) &
422            & .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
423            & CALL ctl_warn( ' Validity time of assimilation increments is ', &
424            &                ' outside the assimilation interval' )
425
426         IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
427            & CALL ctl_warn( ' Validity time of assimilation increments does ', &
428            &                ' not agree with Direct Initialization time' )
429
430         IF ( ln_trainc ) THEN   
431            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
432            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
433            ! Apply the masks
434            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
435            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
436            ! Set missing increments to 0.0 rather than 1e+20
437            ! to allow for differences in masks
438            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
439            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
440         ENDIF
441
442         IF ( ln_dyninc ) THEN   
443            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )             
444            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )             
445            ! Apply the masks
446            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
447            v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
448            ! Set missing increments to 0.0 rather than 1e+20
449            ! to allow for differences in masks
450            WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
451            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
452         ENDIF
453       
454         IF ( ln_sshinc ) THEN
455            CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
456            ! Apply the masks
457            ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
458            ! Set missing increments to 0.0 rather than 1e+20
459            ! to allow for differences in masks
460            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
461         ENDIF
462
463         IF ( ln_seaiceinc ) THEN
464            CALL iom_get( inum, jpdom_autoglo, 'bckinseaice', seaice_bkginc, 1 )
465            ! Apply the masks
466            seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
467            ! Set missing increments to 0.0 rather than 1e+20
468            ! to allow for differences in masks
469            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
470         ENDIF
471
472         IF ( lk_bgcinc ) THEN
473            CALL asm_bgc_init_incs( inum )
474         ENDIF
475
476         CALL iom_close( inum )
477 
478      ENDIF
479
480      !-----------------------------------------------------------------------
481      ! Apply divergence damping filter
482      !-----------------------------------------------------------------------
483
484      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN
485
486         CALL wrk_alloc(jpi,jpj,hdiv) 
487
488         DO  jt = 1, nn_divdmp
489
490            DO jk = 1, jpkm1
491
492               hdiv(:,:) = 0._wp
493
494               DO jj = 2, jpjm1
495                  DO ji = fs_2, fs_jpim1   ! vector opt.
496                     hdiv(ji,jj) =   &
497                        (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * u_bkginc(ji  ,jj  ,jk)     &
498                         - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * u_bkginc(ji-1,jj  ,jk)     &
499                         + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * v_bkginc(ji  ,jj  ,jk)     &
500                         - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * v_bkginc(ji  ,jj-1,jk)  )  &
501                         / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
502                  END DO
503               END DO
504
505               CALL lbc_lnk( hdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
506
507               DO jj = 2, jpjm1
508                  DO ji = fs_2, fs_jpim1   ! vector opt.
509                     u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj)   &
510                                                                        - e1t(ji  ,jj)*e2t(ji  ,jj) * hdiv(ji  ,jj) ) &
511                                                                      / e1u(ji,jj) * umask(ji,jj,jk) 
512                     v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2_wp * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1)   &
513                                                                        - e1t(ji,jj  )*e2t(ji,jj  ) * hdiv(ji,jj  ) ) &
514                                                                      / e2v(ji,jj) * vmask(ji,jj,jk) 
515                  END DO
516               END DO
517
518            END DO
519
520         END DO
521
522         CALL wrk_dealloc(jpi,jpj,hdiv) 
523
524      ENDIF
525
526
527
528      !-----------------------------------------------------------------------
529      ! Allocate and initialize the background state arrays
530      !-----------------------------------------------------------------------
531
532      IF ( ln_asmdin ) THEN
533
534         ALLOCATE( t_bkg(jpi,jpj,jpk) )
535         ALLOCATE( s_bkg(jpi,jpj,jpk) )
536         ALLOCATE( u_bkg(jpi,jpj,jpk) )
537         ALLOCATE( v_bkg(jpi,jpj,jpk) )
538         ALLOCATE( ssh_bkg(jpi,jpj)   )
539
540         t_bkg(:,:,:) = 0.0
541         s_bkg(:,:,:) = 0.0
542         u_bkg(:,:,:) = 0.0
543         v_bkg(:,:,:) = 0.0
544         ssh_bkg(:,:) = 0.0
545
546         !--------------------------------------------------------------------
547         ! Read from file the background state at analysis time
548         !--------------------------------------------------------------------
549
550         CALL iom_open( c_asmdin, inum )
551
552         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
553       
554         IF(lwp) THEN
555            WRITE(numout,*) 
556            WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
557               &  NINT( zdate_bkg )
558            WRITE(numout,*) '~~~~~~~~~~~~'
559         ENDIF
560
561         IF ( NINT( zdate_bkg ) /= iitdin_date ) &
562            & CALL ctl_warn( ' Validity time of assimilation background state does', &
563            &                ' not agree with Direct Initialization time' )
564
565         IF ( ln_trainc ) THEN   
566            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
567            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
568            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
569            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
570         ENDIF
571
572         IF ( ln_dyninc ) THEN   
573            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
574            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
575            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
576            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
577         ENDIF
578       
579         IF ( ln_sshinc ) THEN
580            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
581            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
582         ENDIF
583
584         CALL iom_close( inum )
585
586      ENDIF
587         
588      IF ( lk_bgcinc ) THEN
589         CALL asm_bgc_init_bkg
590      ENDIF
591      !
592   END SUBROUTINE asm_inc_init
593
594
595   SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
596      !!----------------------------------------------------------------------
597      !!                    ***  ROUTINE calc_date  ***
598      !!         
599      !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
600      !!
601      !! ** Method  : Compute the calendar date YYYYMMDD at a given time step.
602      !!
603      !! ** Action  :
604      !!----------------------------------------------------------------------
605      INTEGER, INTENT(IN) :: kit000  ! Initial time step
606      INTEGER, INTENT(IN) :: kt      ! Current time step referenced to kit000
607      INTEGER, INTENT(IN) :: kdate0  ! Initial date
608      INTEGER, INTENT(OUT) :: kdate  ! Current date reference to kdate0
609      !
610      INTEGER :: iyea0    ! Initial year
611      INTEGER :: imon0    ! Initial month
612      INTEGER :: iday0    ! Initial day
613      INTEGER :: iyea     ! Current year
614      INTEGER :: imon     ! Current month
615      INTEGER :: iday     ! Current day
616      INTEGER :: idaystp  ! Number of days between initial and current date
617      INTEGER :: idaycnt  ! Day counter
618
619      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
620
621      !-----------------------------------------------------------------------
622      ! Compute the calendar date YYYYMMDD
623      !-----------------------------------------------------------------------
624
625      ! Initial date
626      iyea0 =   kdate0 / 10000
627      imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
628      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 ) 
629
630      ! Check that kt >= kit000 - 1
631      IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
632
633      ! If kt = kit000 - 1 then set the date to the restart date
634      IF ( kt == kit000 - 1 ) THEN
635
636         kdate = ndastp
637         RETURN
638
639      ENDIF
640
641      ! Compute the number of days from the initial date
642      idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
643   
644      iday    = iday0
645      imon    = imon0
646      iyea    = iyea0
647      idaycnt = 0
648
649      CALL calc_month_len( iyea, imonth_len )
650
651      DO WHILE ( idaycnt < idaystp )
652         iday = iday + 1
653         IF ( iday > imonth_len(imon) )  THEN
654            iday = 1
655            imon = imon + 1
656         ENDIF
657         IF ( imon > 12 ) THEN
658            imon = 1
659            iyea = iyea + 1
660            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
661         ENDIF                 
662         idaycnt = idaycnt + 1
663      END DO
664      !
665      kdate = iyea * 10000 + imon * 100 + iday
666      !
667   END SUBROUTINE
668
669
670   SUBROUTINE calc_month_len( iyear, imonth_len )
671      !!----------------------------------------------------------------------
672      !!                    ***  ROUTINE calc_month_len  ***
673      !!         
674      !! ** Purpose : Compute the number of days in a months given a year.
675      !!
676      !! ** Method  :
677      !!----------------------------------------------------------------------
678      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
679      INTEGER :: iyear         !: year
680      !!----------------------------------------------------------------------
681      !
682      ! length of the month of the current year (from nleapy, read in namelist)
683      IF ( nleapy < 2 ) THEN
684         imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
685         IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years
686            IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
687               imonth_len(2) = 29
688            ENDIF
689         ENDIF
690      ELSE
691         imonth_len(:) = nleapy   ! all months with nleapy days per year
692      ENDIF
693      !
694   END SUBROUTINE
695
696
697   SUBROUTINE tra_asm_inc( kt )
698      !!----------------------------------------------------------------------
699      !!                    ***  ROUTINE tra_asm_inc  ***
700      !!         
701      !! ** Purpose : Apply the tracer (T and S) assimilation increments
702      !!
703      !! ** Method  : Direct initialization or Incremental Analysis Updating
704      !!
705      !! ** Action  :
706      !!----------------------------------------------------------------------
707      INTEGER, INTENT(IN) :: kt               ! Current time step
708      !
709      INTEGER :: ji,jj,jk
710      INTEGER :: it
711      REAL(wp) :: zincwgt  ! IAU weight for current time step
712      REAL (wp), DIMENSION(jpi,jpj,jpk) :: fzptnz ! 3d freezing point values
713      !!----------------------------------------------------------------------
714
715      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
716      ! used to prevent the applied increments taking the temperature below the local freezing point
717
718      DO jk = 1, jpkm1
719        CALL eos_fzp( tsn(:,:,jk,jp_sal), fzptnz(:,:,jk), fsdept(:,:,jk) )
720      END DO
721
722      IF ( ln_asmiau ) THEN
723
724         !--------------------------------------------------------------------
725         ! Incremental Analysis Updating
726         !--------------------------------------------------------------------
727
728         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
729
730            it = kt - nit000 + 1
731            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
732
733            IF(lwp) THEN
734               WRITE(numout,*) 
735               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
736               WRITE(numout,*) '~~~~~~~~~~~~'
737            ENDIF
738
739            ! Update the tracer tendencies
740            DO jk = 1, jpkm1
741               IF (ln_temnofreeze) THEN
742                  ! Do not apply negative increments if the temperature will fall below freezing
743                  WHERE(t_bkginc(:,:,jk) > 0.0_wp .OR. &
744                     &   tsn(:,:,jk,jp_tem) + tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * wgtiau(it) > fzptnz(:,:,jk) ) 
745                     tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
746                  END WHERE
747               ELSE
748                  tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
749               ENDIF
750               IF (ln_salfix) THEN
751                  ! Do not apply negative increments if the salinity will fall below a specified
752                  ! minimum value salfixmin
753                  WHERE(s_bkginc(:,:,jk) > 0.0_wp .OR. &
754                     &   tsn(:,:,jk,jp_sal) + tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * wgtiau(it) > salfixmin ) 
755                     tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
756                  END WHERE
757               ELSE
758                  tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
759               ENDIF
760            END DO
761
762         ENDIF
763
764         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
765            DEALLOCATE( t_bkginc )
766            DEALLOCATE( s_bkginc )
767         ENDIF
768
769
770      ELSEIF ( ln_asmdin ) THEN
771
772         !--------------------------------------------------------------------
773         ! Direct Initialization
774         !--------------------------------------------------------------------
775           
776         IF ( kt == nitdin_r ) THEN
777
778            neuler = 0  ! Force Euler forward step
779
780            ! Initialize the now fields with the background + increment
781            IF (ln_temnofreeze) THEN
782               ! Do not apply negative increments if the temperature will fall below freezing
783               WHERE( t_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_tem) + t_bkginc(:,:,:) > fzptnz(:,:,:) ) 
784                  tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
785               END WHERE
786            ELSE
787               tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
788            ENDIF
789            IF (ln_salfix) THEN
790               ! Do not apply negative increments if the salinity will fall below a specified
791               ! minimum value salfixmin
792               WHERE( s_bkginc(:,:,:) > 0.0_wp .OR. tsn(:,:,:,jp_sal) + s_bkginc(:,:,:) > salfixmin ) 
793                  tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
794               END WHERE
795            ELSE
796               tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
797            ENDIF
798
799            tsb(:,:,:,:) = tsn(:,:,:,:)                 ! Update before fields
800
801            CALL eos( tsb, rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
802!!gm  fabien
803!            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
804!!gm
805
806
807            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)      &
808               &  CALL zps_hde    ( kt, jpts, tsb, gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
809               &                              rhd, gru , grv          )  ! of t, s, rd at the last ocean level
810            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)      &
811               &  CALL zps_hde_isf( nit000, jpts, tsb, gtsu, gtsv,    &    ! Partial steps for top cell (ISF)
812               &                                  rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
813               &                           gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
814
815#if defined key_zdfkpp
816            CALL eos( tsn, rhd, fsdept_n(:,:,:) )                      ! Compute rhd
817!!gm fabien            CALL eos( tsn, rhd )                      ! Compute rhd
818#endif
819
820            DEALLOCATE( t_bkginc )
821            DEALLOCATE( s_bkginc )
822            DEALLOCATE( t_bkg    )
823            DEALLOCATE( s_bkg    )
824         ENDIF
825         
826      ENDIF
827      ! Perhaps the following call should be in step
828      IF   ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
829      !
830   END SUBROUTINE tra_asm_inc
831
832
833   SUBROUTINE dyn_asm_inc( kt )
834      !!----------------------------------------------------------------------
835      !!                    ***  ROUTINE dyn_asm_inc  ***
836      !!         
837      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
838      !!
839      !! ** Method  : Direct initialization or Incremental Analysis Updating.
840      !!
841      !! ** Action  :
842      !!----------------------------------------------------------------------
843      INTEGER, INTENT(IN) :: kt   ! Current time step
844      !
845      INTEGER :: jk
846      INTEGER :: it
847      REAL(wp) :: zincwgt  ! IAU weight for current time step
848      !!----------------------------------------------------------------------
849
850      IF ( ln_asmiau ) THEN
851
852         !--------------------------------------------------------------------
853         ! Incremental Analysis Updating
854         !--------------------------------------------------------------------
855
856         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
857
858            it = kt - nit000 + 1
859            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
860
861            IF(lwp) THEN
862               WRITE(numout,*) 
863               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
864                  &  kt,' with IAU weight = ', wgtiau(it)
865               WRITE(numout,*) '~~~~~~~~~~~~'
866            ENDIF
867
868            ! Update the dynamic tendencies
869            DO jk = 1, jpkm1
870               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
871               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
872            END DO
873           
874            IF ( kt == nitiaufin_r ) THEN
875               DEALLOCATE( u_bkginc )
876               DEALLOCATE( v_bkginc )
877            ENDIF
878
879         ENDIF
880
881      ELSEIF ( ln_asmdin ) THEN 
882
883         !--------------------------------------------------------------------
884         ! Direct Initialization
885         !--------------------------------------------------------------------
886         
887         IF ( kt == nitdin_r ) THEN
888
889            neuler = 0                    ! Force Euler forward step
890
891            ! Initialize the now fields with the background + increment
892            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
893            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
894
895            ub(:,:,:) = un(:,:,:)         ! Update before fields
896            vb(:,:,:) = vn(:,:,:)
897 
898            DEALLOCATE( u_bkg    )
899            DEALLOCATE( v_bkg    )
900            DEALLOCATE( u_bkginc )
901            DEALLOCATE( v_bkginc )
902         ENDIF
903         !
904      ENDIF
905      !
906   END SUBROUTINE dyn_asm_inc
907
908
909   SUBROUTINE ssh_asm_inc( kt )
910      !!----------------------------------------------------------------------
911      !!                    ***  ROUTINE ssh_asm_inc  ***
912      !!         
913      !! ** Purpose : Apply the sea surface height assimilation increment.
914      !!
915      !! ** Method  : Direct initialization or Incremental Analysis Updating.
916      !!
917      !! ** Action  :
918      !!----------------------------------------------------------------------
919      INTEGER, INTENT(IN) :: kt   ! Current time step
920      !
921      INTEGER :: it
922      INTEGER :: jk
923      REAL(wp) :: zincwgt  ! IAU weight for current time step
924      !!----------------------------------------------------------------------
925
926      IF ( ln_asmiau ) THEN
927
928         !--------------------------------------------------------------------
929         ! Incremental Analysis Updating
930         !--------------------------------------------------------------------
931
932         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
933
934            it = kt - nit000 + 1
935            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
936
937            IF(lwp) THEN
938               WRITE(numout,*) 
939               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
940                  &  kt,' with IAU weight = ', wgtiau(it)
941               WRITE(numout,*) '~~~~~~~~~~~~'
942            ENDIF
943
944            ! Save the tendency associated with the IAU weighted SSH increment
945            ! (applied in dynspg.*)
946#if defined key_asminc
947            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
948#endif
949            IF ( kt == nitiaufin_r ) THEN
950               DEALLOCATE( ssh_bkginc )
951            ENDIF
952
953         ELSE
954#if defined key_asminc
955            ssh_iau(:,:) = 0.0
956#endif
957         ENDIF
958
959      ELSEIF ( ln_asmdin ) THEN
960
961         !--------------------------------------------------------------------
962         ! Direct Initialization
963         !--------------------------------------------------------------------
964
965         IF ( kt == nitdin_r ) THEN
966
967            neuler = 0                    ! Force Euler forward step
968
969            ! Initialize the now fields the background + increment
970            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 
971
972            ! Update before fields
973            sshb(:,:) = sshn(:,:)         
974
975            IF( lk_vvl ) THEN
976               DO jk = 1, jpk
977                  fse3t_b(:,:,jk) = fse3t_n(:,:,jk)
978               END DO
979            ENDIF
980
981            DEALLOCATE( ssh_bkg    )
982            DEALLOCATE( ssh_bkginc )
983
984         ENDIF
985         !
986      ENDIF
987      !
988   END SUBROUTINE ssh_asm_inc
989
990
991   SUBROUTINE seaice_asm_inc( kt, kindic )
992      !!----------------------------------------------------------------------
993      !!                    ***  ROUTINE seaice_asm_inc  ***
994      !!         
995      !! ** Purpose : Apply the sea ice assimilation increment.
996      !!
997      !! ** Method  : Direct initialization or Incremental Analysis Updating.
998      !!
999      !! ** Action  :
1000      !!
1001      !!----------------------------------------------------------------------
1002      IMPLICIT NONE
1003      !
1004      INTEGER, INTENT(in)           ::   kt   ! Current time step
1005      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
1006      !
1007      INTEGER  ::   it
1008      REAL(wp) ::   zincwgt   ! IAU weight for current time step
1009#if defined key_lim2
1010      REAL(wp), DIMENSION(jpi,jpj) ::   zofrld, zohicif, zseaicendg, zhicifinc  ! LIM
1011      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
1012#endif
1013      !!----------------------------------------------------------------------
1014
1015      IF ( ln_asmiau ) THEN
1016
1017         !--------------------------------------------------------------------
1018         ! Incremental Analysis Updating
1019         !--------------------------------------------------------------------
1020
1021         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
1022
1023            it = kt - nit000 + 1
1024            zincwgt = wgtiau(it)      ! IAU weight for the current time step
1025            ! note this is not a tendency so should not be divided by rdt (as with the tracer and other increments)
1026
1027            IF(lwp) THEN
1028               WRITE(numout,*) 
1029               WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', &
1030                  &  kt,' with IAU weight = ', wgtiau(it)
1031               WRITE(numout,*) '~~~~~~~~~~~~'
1032            ENDIF
1033
1034            ! Sea-ice : LIM-3 case (to add)
1035
1036#if defined key_lim2
1037            ! Sea-ice : LIM-2 case
1038            zofrld (:,:) = frld(:,:)
1039            zohicif(:,:) = hicif(:,:)
1040            !
1041            frld  = MIN( MAX( frld (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
1042            pfrld = MIN( MAX( pfrld(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
1043            fr_i(:,:) = 1.0_wp - frld(:,:)        ! adjust ice fraction
1044            !
1045            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)   ! find out actual sea ice nudge applied
1046            !
1047            ! Nudge sea ice depth to bring it up to a required minimum depth
1048            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
1049               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt   
1050            ELSEWHERE
1051               zhicifinc(:,:) = 0.0_wp
1052            END WHERE
1053            !
1054            ! nudge ice depth
1055            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1056            phicif(:,:) = phicif(:,:) + zhicifinc(:,:)       
1057            !
1058            ! seaice salinity balancing (to add)
1059#endif
1060
1061#if defined key_cice && defined key_asminc
1062            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1063            ndaice_da(:,:) = seaice_bkginc(:,:) * zincwgt / rdt
1064#endif
1065
1066            IF ( kt == nitiaufin_r ) THEN
1067               DEALLOCATE( seaice_bkginc )
1068            ENDIF
1069
1070         ELSE
1071
1072#if defined key_cice && defined key_asminc
1073            ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1074            ndaice_da(:,:) = 0.0_wp
1075#endif
1076
1077         ENDIF
1078
1079      ELSEIF ( ln_asmdin ) THEN
1080
1081         !--------------------------------------------------------------------
1082         ! Direct Initialization
1083         !--------------------------------------------------------------------
1084
1085         IF ( kt == nitdin_r ) THEN
1086
1087            neuler = 0                    ! Force Euler forward step
1088
1089            ! Sea-ice : LIM-3 case (to add)
1090
1091#if defined key_lim2
1092            ! Sea-ice : LIM-2 case.
1093            zofrld(:,:)=frld(:,:)
1094            zohicif(:,:)=hicif(:,:)
1095            !
1096            ! Initialize the now fields the background + increment
1097            frld (:,:) = MIN( MAX( frld(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
1098            pfrld(:,:) = frld(:,:) 
1099            fr_i (:,:) = 1.0_wp - frld(:,:)                ! adjust ice fraction
1100            zseaicendg(:,:) = zofrld(:,:) - frld(:,:)      ! find out actual sea ice nudge applied
1101            !
1102            ! Nudge sea ice depth to bring it up to a required minimum depth
1103            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hicif(:,:) < zhicifmin ) 
1104               zhicifinc(:,:) = (zhicifmin - hicif(:,:)) * zincwgt   
1105            ELSEWHERE
1106               zhicifinc(:,:) = 0.0_wp
1107            END WHERE
1108            !
1109            ! nudge ice depth
1110            hicif (:,:) = hicif (:,:) + zhicifinc(:,:)
1111            phicif(:,:) = phicif(:,:)       
1112            !
1113            ! seaice salinity balancing (to add)
1114#endif
1115 
1116#if defined key_cice && defined key_asminc
1117            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
1118           ndaice_da(:,:) = seaice_bkginc(:,:) / rdt
1119#endif
1120           IF ( .NOT. PRESENT(kindic) ) THEN
1121              DEALLOCATE( seaice_bkginc )
1122           END IF
1123
1124         ELSE
1125
1126#if defined key_cice && defined key_asminc
1127            ! Sea-ice : CICE case. Zero ice increment tendency into CICE
1128            ndaice_da(:,:) = 0.0_wp
1129#endif
1130         
1131         ENDIF
1132
1133!#if defined defined key_lim2 || defined key_cice
1134!
1135!            IF (ln_seaicebal ) THEN       
1136!             !! balancing salinity increments
1137!             !! simple case from limflx.F90 (doesn't include a mass flux)
1138!             !! assumption is that as ice concentration is reduced or increased
1139!             !! the snow and ice depths remain constant
1140!             !! note that snow is being created where ice concentration is being increased
1141!             !! - could be more sophisticated and
1142!             !! not do this (but would need to alter h_snow)
1143!
1144!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
1145!
1146!             DO jj = 1, jpj
1147!               DO ji = 1, jpi
1148!           ! calculate change in ice and snow mass per unit area
1149!           ! positive values imply adding salt to the ocean (results from ice formation)
1150!           ! fwf : ice formation and melting
1151!
1152!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rdt
1153!
1154!           ! change salinity down to mixed layer depth
1155!                 mld=hmld_kara(ji,jj)
1156!
1157!           ! prevent small mld
1158!           ! less than 10m can cause salinity instability
1159!                 IF (mld < 10) mld=10
1160!
1161!           ! set to bottom of a level
1162!                 DO jk = jpk-1, 2, -1
1163!                   IF ((mld > gdepw(ji,jj,jk)) .and. (mld < gdepw(ji,jj,jk+1))) THEN
1164!                     mld=gdepw(ji,jj,jk+1)
1165!                     jkmax=jk
1166!                   ENDIF
1167!                 ENDDO
1168!
1169!            ! avoid applying salinity balancing in shallow water or on land
1170!            !
1171!
1172!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
1173!
1174!                 dsal_ocn=0.0_wp
1175!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
1176!
1177!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
1178!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
1179!
1180!           ! put increments in for levels in the mixed layer
1181!           ! but prevent salinity below a threshold value
1182!
1183!                   DO jk = 1, jkmax             
1184!
1185!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1186!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1187!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1188!                     ENDIF
1189!
1190!                   ENDDO
1191!
1192!      !            !  salt exchanges at the ice/ocean interface
1193!      !            zpmess         = zfons / rdt_ice    ! rdt_ice is ice timestep
1194!      !
1195!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1196!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
1197!      !!               
1198!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
1199!      !!                                                     ! E-P (kg m-2 s-2)
1200!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
1201!               ENDDO !ji
1202!             ENDDO !jj!
1203!
1204!            ENDIF !ln_seaicebal
1205!
1206!#endif
1207
1208      ENDIF
1209
1210   END SUBROUTINE seaice_asm_inc
1211
1212
1213   SUBROUTINE bgc_asm_inc( kt )
1214      !!----------------------------------------------------------------------
1215      !!                    ***  ROUTINE bgc_asm_inc  ***
1216      !!         
1217      !! ** Purpose : Apply the biogeochemistry assimilation increments
1218      !!
1219      !! ** Method  : Call relevant routines in asmbgc
1220      !!
1221      !! ** Action  : Call relevant routines in asmbgc
1222      !!
1223      !!----------------------------------------------------------------------
1224      !!
1225      INTEGER, INTENT(in   ) :: kt        ! Current time step
1226      !
1227      INTEGER                :: icycper   ! Dimension of wgtiau
1228      !!
1229      !!----------------------------------------------------------------------
1230     
1231      icycper = SIZE( wgtiau )
1232     
1233      ! Ocean colour variables first
1234      IF ( ln_slchltotinc .OR. ln_slchldiainc .OR. ln_slchlnoninc .OR. &
1235         & ln_schltotinc  .OR. ln_slphytotinc .OR. ln_slphydiainc .OR. &
1236         & ln_slphynoninc ) THEN
1237         CALL phyto2d_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau )
1238      ENDIF
1239     
1240      ! Surface pCO2/fCO2 next
1241      IF ( ln_sfco2inc .OR. ln_spco2inc ) THEN
1242         CALL pco2_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, &
1243            &               ln_trainc, t_bkginc, s_bkginc )
1244      ENDIF
1245     
1246      ! Profile pH next
1247      IF ( ln_pphinc ) THEN
1248         CALL ph_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau, &
1249            &             ln_trainc, t_bkginc, s_bkginc )
1250      ENDIF
1251     
1252      ! Then chlorophyll profiles
1253      IF ( ln_plchltotinc .OR. ln_pchltotinc ) THEN
1254         CALL phyto3d_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau )
1255      ENDIF
1256     
1257      ! Remaining bgc profile variables
1258      IF ( ln_pno3inc .OR. ln_psi4inc .OR. ln_pdicinc .OR. &
1259         & ln_palkinc .OR. ln_po2inc ) THEN
1260         CALL bgc3d_asm_inc( kt, ln_asmdin, ln_asmiau, icycper, wgtiau )
1261      ENDIF
1262
1263   END SUBROUTINE bgc_asm_inc
1264   
1265   !!======================================================================
1266END MODULE asminc
Note: See TracBrowser for help on using the repository browser.