source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_sit/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 12595

Last change on this file since 12595 was 12595, checked in by dcarneir, 8 months ago

Removing commented code in subroutine sit_asm_inc

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