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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
asminc.F90 in NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ASM – NEMO

source: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/ASM/asminc.F90 @ 14002

Last change on this file since 14002 was 14002, checked in by emanuelaclementi, 3 years ago

phasing with trunk rev14001 - ticket #2152 #2339

  • Property svn:keywords set to Id
File size: 51.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   !!                 ! 2014-09  (D. Lea)  Local calc_date removed use routine from OBS
14   !!                 ! 2015-11  (D. Lea)  Handle non-zero initial time of day
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!   asm_inc_init   : Initialize the increment arrays and IAU weights
19   !!   tra_asm_inc    : Apply the tracer (T and S) increments
20   !!   dyn_asm_inc    : Apply the dynamic (u and v) increments
21   !!   ssh_asm_inc    : Apply the SSH increment
22   !!   ssh_asm_div    : Apply divergence associated with SSH increment
23   !!   seaice_asm_inc : Apply the seaice increment
24   !!----------------------------------------------------------------------
25   USE oce             ! Dynamics and active tracers defined in memory
26   USE par_oce         ! Ocean space and time domain variables
27   USE dom_oce         ! Ocean space and time domain
28   USE domain, ONLY : dom_tile
29   USE domvvl          ! domain: variable volume level
30   USE ldfdyn          ! lateral diffusion: eddy viscosity coefficients
31   USE eosbn2          ! Equation of state - in situ and potential density
32   USE zpshde          ! Partial step : Horizontal Derivative
33   USE asmpar          ! Parameters for the assmilation interface
34   USE asmbkg          !
35   USE c1d             ! 1D initialization
36   USE sbc_oce         ! Surface boundary condition variables.
37   USE diaobs   , ONLY : calc_date     ! Compute the calendar date on a given step
38#if defined key_si3
39   USE ice      , ONLY : hm_i, at_i, at_i_b
40#endif
41   !
42   USE in_out_manager  ! I/O manager
43   USE iom             ! Library to read input files
44   USE lib_mpp         ! MPP library
45
46   IMPLICIT NONE
47   PRIVATE
48   
49   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
50   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
51   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
52   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
53   PUBLIC   ssh_asm_div    !: Apply the SSH divergence
54   PUBLIC   seaice_asm_inc !: Apply the seaice increment
55
56#if defined key_asminc
57    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
58#else
59    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
60#endif
61   LOGICAL, PUBLIC :: ln_bkgwri     !: No output of the background state fields
62   LOGICAL, PUBLIC :: ln_asmiau     !: No applying forcing with an assimilation increment
63   LOGICAL, PUBLIC :: ln_asmdin     !: No direct initialization
64   LOGICAL, PUBLIC :: ln_trainc     !: No tracer (T and S) assimilation increments
65   LOGICAL, PUBLIC :: ln_dyninc     !: No dynamics (u and v) assimilation increments
66   LOGICAL, PUBLIC :: ln_sshinc     !: No sea surface height assimilation increment
67   LOGICAL, PUBLIC :: ln_seaiceinc  !: No sea ice concentration increment
68   LOGICAL, PUBLIC :: ln_salfix     !: Apply minimum salinity check
69   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing
70   INTEGER, PUBLIC :: nn_divdmp     !: Apply divergence damping filter nn_divdmp times
71
72   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
73   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
74   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
75   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
76   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
77#if defined key_asminc
78   REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE ::   ssh_iau              !: IAU-weighted sea surface height increment
79#endif
80   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
81   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
82   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
83   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
84   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
85   !
86   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
87   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
88   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
89
90   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
91   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc
92#if defined key_cice && defined key_asminc
93   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ndaice_da             ! ice increment tendency into CICE
94#endif
95
96   !! * Substitutions
97#  include "do_loop_substitute.h90"
98#  include "domzgr_substitute.h90"
99   !!----------------------------------------------------------------------
100   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
101   !! $Id$
102   !! Software governed by the CeCILL license (see ./LICENSE)
103   !!----------------------------------------------------------------------
104CONTAINS
105
106   SUBROUTINE asm_inc_init( Kbb, Kmm, Krhs )
107      !!----------------------------------------------------------------------
108      !!                    ***  ROUTINE asm_inc_init  ***
109      !!         
110      !! ** Purpose : Initialize the assimilation increment and IAU weights.
111      !!
112      !! ** Method  : Initialize the assimilation increment and IAU weights.
113      !!
114      !! ** Action  :
115      !!----------------------------------------------------------------------
116      INTEGER, INTENT(in) ::  Kbb, Kmm, Krhs  ! time level indices
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      REAL(KIND=dp) :: ditend_date     ! Date YYYYMMDD.HHMMSS of final time step
124      REAL(KIND=dp) :: ditbkg_date     ! Date YYYYMMDD.HHMMSS of background time step for Jb term
125      REAL(KIND=dp) :: ditdin_date     ! Date YYYYMMDD.HHMMSS of background time step for DI
126      REAL(KIND=dp) :: ditiaustr_date  ! Date YYYYMMDD.HHMMSS of IAU interval start time step
127      REAL(KIND=dp) :: ditiaufin_date  ! Date YYYYMMDD.HHMMSS 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), ALLOCATABLE, DIMENSION(:,:) ::   zhdiv   ! 2D workspace
137      !!
138      NAMELIST/nam_asminc/ ln_bkgwri,                                      &
139         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
140         &                 ln_asmdin, ln_asmiau,                           &
141         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
142         &                 ln_salfix, salfixmin, nn_divdmp
143      !!----------------------------------------------------------------------
144
145      !-----------------------------------------------------------------------
146      ! Read Namelist nam_asminc : assimilation increment interface
147      !-----------------------------------------------------------------------
148      ln_seaiceinc   = .FALSE.
149      ln_temnofreeze = .FALSE.
150
151      READ  ( numnam_ref, nam_asminc, IOSTAT = ios, ERR = 901)
152901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_asminc in reference namelist' )
153      READ  ( numnam_cfg, nam_asminc, IOSTAT = ios, ERR = 902 )
154902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_asminc in configuration namelist' )
155      IF(lwm) WRITE ( numond, nam_asminc )
156
157      ! Control print
158      IF(lwp) THEN
159         WRITE(numout,*)
160         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
161         WRITE(numout,*) '~~~~~~~~~~~~'
162         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
163         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
164         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
165         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
166         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
167         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
168         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc
169         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
170         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
171         WRITE(numout,*) '      Timestep of background for DI in [0,nitend-nit000-1]     nitdin    = ', nitdin
172         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
173         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
174         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
175         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
176         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
177      ENDIF
178
179      nitbkg_r    = nitbkg    + nit000 - 1            ! Background time referenced to nit000
180      nitdin_r    = nitdin    + nit000 - 1            ! Background time for DI referenced to nit000
181      nitiaustr_r = nitiaustr + nit000 - 1            ! Start of IAU interval referenced to nit000
182      nitiaufin_r = nitiaufin + nit000 - 1            ! End of IAU interval referenced to nit000
183
184      iiauper     = nitiaufin_r - nitiaustr_r + 1     ! IAU interval length
185      icycper     = nitend      - nit000      + 1     ! Cycle interval length
186
187      CALL calc_date( nitend     , ditend_date    )   ! Date of final time step
188      CALL calc_date( nitbkg_r   , ditbkg_date    )   ! Background time for Jb referenced to ndate0
189      CALL calc_date( nitdin_r   , ditdin_date    )   ! Background time for DI referenced to ndate0
190      CALL calc_date( nitiaustr_r, ditiaustr_date )   ! IAU start time referenced to ndate0
191      CALL calc_date( nitiaufin_r, ditiaufin_date )   ! IAU end time referenced to ndate0
192
193      IF(lwp) THEN
194         WRITE(numout,*)
195         WRITE(numout,*) '   Time steps referenced to current cycle:'
196         WRITE(numout,*) '       iitrst      = ', nit000 - 1
197         WRITE(numout,*) '       nit000      = ', nit000
198         WRITE(numout,*) '       nitend      = ', nitend
199         WRITE(numout,*) '       nitbkg_r    = ', nitbkg_r
200         WRITE(numout,*) '       nitdin_r    = ', nitdin_r
201         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
202         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
203         WRITE(numout,*)
204         WRITE(numout,*) '   Dates referenced to current cycle:'
205         WRITE(numout,*) '       ndastp         = ', ndastp
206         WRITE(numout,*) '       ndate0         = ', ndate0
207         WRITE(numout,*) '       nn_time0       = ', nn_time0
208         WRITE(numout,*) '       ditend_date    = ', ditend_date
209         WRITE(numout,*) '       ditbkg_date    = ', ditbkg_date
210         WRITE(numout,*) '       ditdin_date    = ', ditdin_date
211         WRITE(numout,*) '       ditiaustr_date = ', ditiaustr_date
212         WRITE(numout,*) '       ditiaufin_date = ', ditiaufin_date
213      ENDIF
214
215
216      IF ( ( ln_asmdin ).AND.( ln_asmiau ) )   &
217         & CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
218         &                ' Choose Direct Initialization OR Incremental Analysis Updating')
219
220      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
221           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) &
222         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', &
223         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
224         &                ' Inconsistent options')
225
226      IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
227         & CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :',  &
228         &                ' Type IAU weighting function is invalid')
229
230      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) &
231         &                     )  &
232         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', &
233         &                ' The assimilation increments are not applied')
234
235      IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
236         & CALL ctl_stop( ' nitiaustr = nitiaufin :',  &
237         &                ' IAU interval is of zero length')
238
239      IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
240         & CALL ctl_stop( ' nitiaustr or nitiaufin :',  &
241         &                ' IAU starting or final time step is outside the cycle interval', &
242         &                 ' Valid range nit000 to nitend')
243
244      IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
245         & CALL ctl_stop( ' nitbkg :',  &
246         &                ' Background time step is outside the cycle interval')
247
248      IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
249         & CALL ctl_stop( ' nitdin :',  &
250         &                ' Background time step for Direct Initialization is outside', &
251         &                ' the cycle interval')
252
253      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
254
255      !--------------------------------------------------------------------
256      ! Initialize the Incremental Analysis Updating weighting function
257      !--------------------------------------------------------------------
258
259      IF( ln_asmiau ) THEN
260         !
261         ALLOCATE( wgtiau( icycper ) )
262         !
263         wgtiau(:) = 0._wp
264         !
265         !                                !---------------------------------------------------------
266         IF( niaufn == 0 ) THEN           ! Constant IAU forcing
267            !                             !---------------------------------------------------------
268            DO jt = 1, iiauper
269               wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
270            END DO
271            !                             !---------------------------------------------------------
272         ELSEIF ( niaufn == 1 ) THEN      ! Linear hat-like, centred in middle of IAU interval
273            !                             !---------------------------------------------------------
274            ! Compute the normalization factor
275            znorm = 0._wp
276            IF( MOD( iiauper, 2 ) == 0 ) THEN   ! Even number of time steps in IAU interval
277               imid = iiauper / 2 
278               DO jt = 1, imid
279                  znorm = znorm + REAL( jt )
280               END DO
281               znorm = 2.0 * znorm
282            ELSE                                ! Odd number of time steps in IAU interval
283               imid = ( iiauper + 1 ) / 2       
284               DO jt = 1, imid - 1
285                  znorm = znorm + REAL( jt )
286               END DO
287               znorm = 2.0 * znorm + REAL( imid )
288            ENDIF
289            znorm = 1.0 / znorm
290            !
291            DO jt = 1, imid - 1
292               wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
293            END DO
294            DO jt = imid, iiauper
295               wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
296            END DO
297            !
298         ENDIF
299
300         ! Test that the integral of the weights over the weighting interval equals 1
301          IF(lwp) THEN
302             WRITE(numout,*)
303             WRITE(numout,*) 'asm_inc_init : IAU weights'
304             WRITE(numout,*) '~~~~~~~~~~~~'
305             WRITE(numout,*) '             time step         IAU  weight'
306             WRITE(numout,*) '             =========     ====================='
307             ztotwgt = 0.0
308             DO jt = 1, icycper
309                ztotwgt = ztotwgt + wgtiau(jt)
310                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
311             END DO   
312             WRITE(numout,*) '         ==================================='
313             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
314             WRITE(numout,*) '         ==================================='
315          ENDIF
316         
317      ENDIF
318
319      !--------------------------------------------------------------------
320      ! Allocate and initialize the increment arrays
321      !--------------------------------------------------------------------
322
323      ALLOCATE( t_bkginc     (jpi,jpj,jpk) )   ;   t_bkginc     (:,:,:) = 0._wp
324      ALLOCATE( s_bkginc     (jpi,jpj,jpk) )   ;   s_bkginc     (:,:,:) = 0._wp
325      ALLOCATE( u_bkginc     (jpi,jpj,jpk) )   ;   u_bkginc     (:,:,:) = 0._wp
326      ALLOCATE( v_bkginc     (jpi,jpj,jpk) )   ;   v_bkginc     (:,:,:) = 0._wp
327      ALLOCATE( ssh_bkginc   (jpi,jpj)     )   ;   ssh_bkginc   (:,:)   = 0._wp
328      ALLOCATE( seaice_bkginc(jpi,jpj)     )   ;   seaice_bkginc(:,:)   = 0._wp
329#if defined key_asminc
330      ALLOCATE( ssh_iau      (jpi,jpj)     )   ;   ssh_iau      (:,:)   = 0._wp
331#endif
332#if defined key_cice && defined key_asminc
333      ALLOCATE( ndaice_da    (jpi,jpj)     )   ;   ndaice_da    (:,:)   = 0._wp
334#endif
335      !
336      IF ( ln_trainc .OR. ln_dyninc .OR.   &       !--------------------------------------
337         & ln_sshinc .OR. ln_seaiceinc   ) THEN    ! Read the increments from file
338         !                                         !--------------------------------------
339         CALL iom_open( c_asminc, inum )
340         !
341         CALL iom_get( inum, 'time'       , zdate_inc   ) 
342         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
343         CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
344         z_inc_dateb = zdate_inc
345         z_inc_datef = zdate_inc
346         !
347         IF(lwp) THEN
348            WRITE(numout,*) 
349            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid between dates ', z_inc_dateb,' and ', z_inc_datef
350            WRITE(numout,*) '~~~~~~~~~~~~'
351         ENDIF
352         !
353         IF ( ( z_inc_dateb < ndastp + nn_time0*0.0001_wp ) .OR.   &
354            & ( z_inc_datef > ditend_date ) ) &
355            &    CALL ctl_warn( ' Validity time of assimilation increments is ', &
356            &                   ' outside the assimilation interval' )
357
358         IF ( ( ln_asmdin ).AND.( zdate_inc /= ditdin_date ) ) &
359            & CALL ctl_warn( ' Validity time of assimilation increments does ', &
360            &                ' not agree with Direct Initialization time' )
361
362         IF ( ln_trainc ) THEN   
363            CALL iom_get( inum, jpdom_auto, 'bckint', t_bkginc, 1 )
364            CALL iom_get( inum, jpdom_auto, 'bckins', s_bkginc, 1 )
365            ! Apply the masks
366            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
367            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
368            ! Set missing increments to 0.0 rather than 1e+20
369            ! to allow for differences in masks
370            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
371            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
372         ENDIF
373
374         IF ( ln_dyninc ) THEN   
375            CALL iom_get( inum, jpdom_auto, 'bckinu', u_bkginc, 1 )             
376            CALL iom_get( inum, jpdom_auto, 'bckinv', v_bkginc, 1 )             
377            ! Apply the masks
378            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
379            v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
380            ! Set missing increments to 0.0 rather than 1e+20
381            ! to allow for differences in masks
382            WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
383            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
384         ENDIF
385       
386         IF ( ln_sshinc ) THEN
387            CALL iom_get( inum, jpdom_auto, 'bckineta', ssh_bkginc, 1 )
388            ! Apply the masks
389            ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
390            ! Set missing increments to 0.0 rather than 1e+20
391            ! to allow for differences in masks
392            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
393         ENDIF
394
395         IF ( ln_seaiceinc ) THEN
396            CALL iom_get( inum, jpdom_auto, 'bckinseaice', seaice_bkginc, 1 )
397            ! Apply the masks
398            seaice_bkginc(:,:) = seaice_bkginc(:,:) * tmask(:,:,1)
399            ! Set missing increments to 0.0 rather than 1e+20
400            ! to allow for differences in masks
401            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0
402         ENDIF
403         !
404         CALL iom_close( inum )
405         !
406      ENDIF
407      !
408      !                                            !--------------------------------------
409      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN    ! Apply divergence damping filter
410         !                                         !--------------------------------------
411         ALLOCATE( zhdiv(jpi,jpj) ) 
412         !
413         DO jt = 1, nn_divdmp
414            !
415            DO jk = 1, jpkm1           ! zhdiv = e1e1 * div
416               zhdiv(:,:) = 0._wp
417               DO_2D( 0, 0, 0, 0 )
418                  zhdiv(ji,jj) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * u_bkginc(ji  ,jj,jk)    &
419                     &            - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * u_bkginc(ji-1,jj,jk)    &
420                     &            + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * v_bkginc(ji,jj  ,jk)    &
421                     &            - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * v_bkginc(ji,jj-1,jk)  ) &
422                     &            / e3t(ji,jj,jk,Kmm)
423               END_2D
424               CALL lbc_lnk( 'asminc', zhdiv, 'T', 1.0_wp )   ! lateral boundary cond. (no sign change)
425               !
426               DO_2D( 0, 0, 0, 0 )
427                  u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk)                         &
428                     &               + 0.2_wp * ( zhdiv(ji+1,jj) - zhdiv(ji  ,jj) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
429                  v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk)                         &
430                     &               + 0.2_wp * ( zhdiv(ji,jj+1) - zhdiv(ji,jj  ) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 
431               END_2D
432            END DO
433            !
434         END DO
435         !
436         DEALLOCATE( zhdiv ) 
437         !
438      ENDIF
439      !
440      !                             !-----------------------------------------------------
441      IF ( ln_asmdin ) THEN         ! Allocate and initialize the background state arrays
442         !                          !-----------------------------------------------------
443         !
444         ALLOCATE( t_bkg  (jpi,jpj,jpk) )   ;   t_bkg  (:,:,:) = 0._wp
445         ALLOCATE( s_bkg  (jpi,jpj,jpk) )   ;   s_bkg  (:,:,:) = 0._wp
446         ALLOCATE( u_bkg  (jpi,jpj,jpk) )   ;   u_bkg  (:,:,:) = 0._wp
447         ALLOCATE( v_bkg  (jpi,jpj,jpk) )   ;   v_bkg  (:,:,:) = 0._wp
448         ALLOCATE( ssh_bkg(jpi,jpj)     )   ;   ssh_bkg(:,:)   = 0._wp
449         !
450         !
451         !--------------------------------------------------------------------
452         ! Read from file the background state at analysis time
453         !--------------------------------------------------------------------
454         !
455         CALL iom_open( c_asmdin, inum )
456         !
457         CALL iom_get( inum, 'rdastp', zdate_bkg ) 
458         !
459         IF(lwp) THEN
460            WRITE(numout,*) 
461            WRITE(numout,*) '   ==>>>  Assimilation background state valid at : ', zdate_bkg
462            WRITE(numout,*)
463         ENDIF
464         !
465         IF ( zdate_bkg /= ditdin_date )   &
466            & CALL ctl_warn( ' Validity time of assimilation background state does', &
467            &                ' not agree with Direct Initialization time' )
468         !
469         IF ( ln_trainc ) THEN   
470            CALL iom_get( inum, jpdom_auto, 'tn', t_bkg )
471            CALL iom_get( inum, jpdom_auto, 'sn', s_bkg )
472            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
473            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
474         ENDIF
475         !
476         IF ( ln_dyninc ) THEN   
477            CALL iom_get( inum, jpdom_auto, 'un', u_bkg, cd_type = 'U', psgn = 1._wp )
478            CALL iom_get( inum, jpdom_auto, 'vn', v_bkg, cd_type = 'V', psgn = 1._wp )
479            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
480            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
481         ENDIF
482         !
483         IF ( ln_sshinc ) THEN
484            CALL iom_get( inum, jpdom_auto, 'sshn', ssh_bkg )
485            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
486         ENDIF
487         !
488         CALL iom_close( inum )
489         !
490      ENDIF
491      !
492      IF(lwp) WRITE(numout,*) '   ==>>>   Euler time step switch is ', l_1st_euler
493      !
494      IF( lk_asminc ) THEN                            !==  data assimilation  ==!
495         IF( ln_bkgwri )   CALL asm_bkg_wri( nit000 - 1, Kmm )      ! Output background fields
496         IF( ln_asmdin ) THEN                                  ! Direct initialization
497            IF( ln_trainc )   CALL tra_asm_inc( nit000 - 1, Kbb, Kmm, ts    , Krhs )      ! Tracers
498            IF( ln_dyninc )   CALL dyn_asm_inc( nit000 - 1, Kbb, Kmm, uu, vv, Krhs )      ! Dynamics
499            IF( ln_sshinc )   CALL ssh_asm_inc( nit000 - 1, Kbb, Kmm )                    ! SSH
500         ENDIF
501      ENDIF
502      !
503   END SUBROUTINE asm_inc_init
504   
505   
506   SUBROUTINE tra_asm_inc( kt, Kbb, Kmm, pts, Krhs )
507      !!----------------------------------------------------------------------
508      !!                    ***  ROUTINE tra_asm_inc  ***
509      !!         
510      !! ** Purpose : Apply the tracer (T and S) assimilation increments
511      !!
512      !! ** Method  : Direct initialization or Incremental Analysis Updating
513      !!
514      !! ** Action  :
515      !!----------------------------------------------------------------------
516      INTEGER                                  , INTENT(in   ) :: kt             ! Current time step
517      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Krhs ! Time level indices
518      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers and RHS of tracer equation
519      !
520      INTEGER  :: ji, jj, jk
521      INTEGER  :: it, itile
522      REAL(wp) :: zincwgt  ! IAU weight for current time step
523      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: fzptnz ! 3d freezing point values
524      !!----------------------------------------------------------------------
525      !
526      ! freezing point calculation taken from oc_fz_pt (but calculated for all depths)
527      ! used to prevent the applied increments taking the temperature below the local freezing point
528      IF( ln_temnofreeze ) THEN
529         DO jk = 1, jpkm1
530           CALL eos_fzp( pts(:,:,jk,jp_sal,Kmm), fzptnz(:,:,jk), gdept(:,:,jk,Kmm) )
531         END DO
532      ENDIF
533         !
534         !                             !--------------------------------------
535      IF ( ln_asmiau ) THEN            ! Incremental Analysis Updating
536         !                             !--------------------------------------
537         !
538         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
539            !
540            it = kt - nit000 + 1
541            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
542            !
543            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
544               IF(lwp) THEN
545                  WRITE(numout,*)
546                  WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
547                  WRITE(numout,*) '~~~~~~~~~~~~'
548               ENDIF
549            ENDIF
550            !
551            ! Update the tracer tendencies
552            DO jk = 1, jpkm1
553               IF (ln_temnofreeze) THEN
554                  ! Do not apply negative increments if the temperature will fall below freezing
555                  WHERE(t_bkginc(A2D(0),jk) > 0.0_wp .OR. &
556                     &   pts(A2D(0),jk,jp_tem,Kmm) + pts(A2D(0),jk,jp_tem,Krhs) + t_bkginc(A2D(0),jk) * wgtiau(it) > fzptnz(:,:,jk) )
557                     pts(A2D(0),jk,jp_tem,Krhs) = pts(A2D(0),jk,jp_tem,Krhs) + t_bkginc(A2D(0),jk) * zincwgt
558                  END WHERE
559               ELSE
560                  DO_2D( 0, 0, 0, 0 )
561                     pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) + t_bkginc(ji,jj,jk) * zincwgt
562                  END_2D
563               ENDIF
564               IF (ln_salfix) THEN
565                  ! Do not apply negative increments if the salinity will fall below a specified
566                  ! minimum value salfixmin
567                  WHERE(s_bkginc(A2D(0),jk) > 0.0_wp .OR. &
568                     &   pts(A2D(0),jk,jp_sal,Kmm) + pts(A2D(0),jk,jp_sal,Krhs) + s_bkginc(A2D(0),jk) * wgtiau(it) > salfixmin )
569                     pts(A2D(0),jk,jp_sal,Krhs) = pts(A2D(0),jk,jp_sal,Krhs) + s_bkginc(A2D(0),jk) * zincwgt
570                  END WHERE
571               ELSE
572                  DO_2D( 0, 0, 0, 0 )
573                     pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) + s_bkginc(ji,jj,jk) * zincwgt
574                  END_2D
575               ENDIF
576            END DO
577            !
578         ENDIF
579         !
580         IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
581            IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
582               DEALLOCATE( t_bkginc )
583               DEALLOCATE( s_bkginc )
584            ENDIF
585         ENDIF
586         !                             !--------------------------------------
587      ELSEIF ( ln_asmdin ) THEN        ! Direct Initialization
588         !                             !--------------------------------------
589         !           
590         IF ( kt == nitdin_r ) THEN
591            !
592            l_1st_euler = .TRUE.  ! Force Euler forward step
593            !
594            ! Initialize the now fields with the background + increment
595            IF (ln_temnofreeze) THEN
596               ! Do not apply negative increments if the temperature will fall below freezing
597               WHERE( t_bkginc(A2D(0),:) > 0.0_wp .OR. pts(A2D(0),:,jp_tem,Kmm) + t_bkginc(A2D(0),:) > fzptnz(:,:,:) )
598                  pts(A2D(0),:,jp_tem,Kmm) = t_bkg(A2D(0),:) + t_bkginc(A2D(0),:)
599               END WHERE
600            ELSE
601               DO_3D( 0, 0, 0, 0, 1, jpk )
602                  pts(ji,jj,jk,jp_tem,Kmm) = t_bkg(ji,jj,jk) + t_bkginc(ji,jj,jk)
603               END_3D
604            ENDIF
605            IF (ln_salfix) THEN
606               ! Do not apply negative increments if the salinity will fall below a specified
607               ! minimum value salfixmin
608               WHERE( s_bkginc(A2D(0),:) > 0.0_wp .OR. pts(A2D(0),:,jp_sal,Kmm) + s_bkginc(A2D(0),:) > salfixmin )
609                  pts(A2D(0),:,jp_sal,Kmm) = s_bkg(A2D(0),:) + s_bkginc(A2D(0),:)
610               END WHERE
611            ELSE
612               DO_3D( 0, 0, 0, 0, 1, jpk )
613                  pts(ji,jj,jk,jp_sal,Kmm) = s_bkg(ji,jj,jk) + s_bkginc(ji,jj,jk)
614               END_3D
615            ENDIF
616
617            DO_3D( 0, 0, 0, 0, 1, jpk )
618               pts(ji,jj,jk,:,Kbb) = pts(ji,jj,jk,:,Kmm)             ! Update before fields
619            END_3D
620
621            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop, gdept_0(:,:,:) )  ! Before potential and in situ densities
622!!gm  fabien
623!            CALL eos( pts(:,:,:,:,Kbb), rhd, rhop )                ! Before potential and in situ densities
624!!gm
625
626            ! TEMP: [tiling] This change not necessary after extra haloes development (lbc_lnk removed from zps_hde*)
627            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only for the full domain
628               itile = ntile
629               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )            ! Use full domain
630
631               IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav)           &
632                  &  CALL zps_hde    ( kt, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv,        &  ! Partial steps: before horizontal gradient
633                  &                              rhd, gru , grv               )  ! of t, s, rd at the last ocean level
634               IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav)                       &
635                  &  CALL zps_hde_isf( nit000, Kmm, jpts, pts(:,:,:,:,Kbb), gtsu, gtsv, gtui, gtvi,    &  ! Partial steps for top cell (ISF)
636                  &                                  rhd, gru , grv , grui, grvi          )  ! of t, s, rd at the last ocean level
637
638               IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = itile )            ! Revert to tile domain
639            ENDIF
640
641            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
642               DEALLOCATE( t_bkginc )
643               DEALLOCATE( s_bkginc )
644               DEALLOCATE( t_bkg    )
645               DEALLOCATE( s_bkg    )
646            ENDIF
647         !
648         ENDIF
649         
650      ENDIF
651      ! Perhaps the following call should be in step
652      IF ( ln_seaiceinc  )   CALL seaice_asm_inc ( kt )   ! apply sea ice concentration increment
653      !
654   END SUBROUTINE tra_asm_inc
655
656
657   SUBROUTINE dyn_asm_inc( kt, Kbb, Kmm, puu, pvv, Krhs )
658      !!----------------------------------------------------------------------
659      !!                    ***  ROUTINE dyn_asm_inc  ***
660      !!         
661      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
662      !!
663      !! ** Method  : Direct initialization or Incremental Analysis Updating.
664      !!
665      !! ** Action  :
666      !!----------------------------------------------------------------------
667      INTEGER                             , INTENT( in )  ::  kt             ! ocean time-step index
668      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs ! ocean time level indices
669      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv       ! ocean velocities and RHS of momentum equation
670      !
671      INTEGER :: jk
672      INTEGER :: it
673      REAL(wp) :: zincwgt  ! IAU weight for current time step
674      !!----------------------------------------------------------------------
675      !
676      !                          !--------------------------------------------
677      IF ( ln_asmiau ) THEN      ! Incremental Analysis Updating
678         !                       !--------------------------------------------
679         !
680         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
681            !
682            it = kt - nit000 + 1
683            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
684            !
685            IF(lwp) THEN
686               WRITE(numout,*) 
687               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
688               WRITE(numout,*) '~~~~~~~~~~~~'
689            ENDIF
690            !
691            ! Update the dynamic tendencies
692            DO jk = 1, jpkm1
693               puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + u_bkginc(:,:,jk) * zincwgt
694               pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + v_bkginc(:,:,jk) * zincwgt
695            END DO
696            !
697            IF ( kt == nitiaufin_r ) THEN
698               DEALLOCATE( u_bkginc )
699               DEALLOCATE( v_bkginc )
700            ENDIF
701            !
702         ENDIF
703         !                          !-----------------------------------------
704      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
705         !                          !-----------------------------------------
706         !         
707         IF ( kt == nitdin_r ) THEN
708            !
709            l_1st_euler = .TRUE.                    ! Force Euler forward step
710            !
711            ! Initialize the now fields with the background + increment
712            puu(:,:,:,Kmm) = u_bkg(:,:,:) + u_bkginc(:,:,:)
713            pvv(:,:,:,Kmm) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
714            !
715            puu(:,:,:,Kbb) = puu(:,:,:,Kmm)         ! Update before fields
716            pvv(:,:,:,Kbb) = pvv(:,:,:,Kmm)
717            !
718            DEALLOCATE( u_bkg    )
719            DEALLOCATE( v_bkg    )
720            DEALLOCATE( u_bkginc )
721            DEALLOCATE( v_bkginc )
722         ENDIF
723         !
724      ENDIF
725      !
726   END SUBROUTINE dyn_asm_inc
727
728
729   SUBROUTINE ssh_asm_inc( kt, Kbb, Kmm )
730      !!----------------------------------------------------------------------
731      !!                    ***  ROUTINE ssh_asm_inc  ***
732      !!         
733      !! ** Purpose : Apply the sea surface height assimilation increment.
734      !!
735      !! ** Method  : Direct initialization or Incremental Analysis Updating.
736      !!
737      !! ** Action  :
738      !!----------------------------------------------------------------------
739      INTEGER, INTENT(IN) :: kt         ! Current time step
740      INTEGER, INTENT(IN) :: Kbb, Kmm   ! Current time step
741      !
742      INTEGER :: it
743      INTEGER :: jk
744      REAL(wp) :: zincwgt  ! IAU weight for current time step
745      !!----------------------------------------------------------------------
746      !
747      !                             !-----------------------------------------
748      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
749         !                          !-----------------------------------------
750         !
751         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
752            !
753            it = kt - nit000 + 1
754            zincwgt = wgtiau(it) / rn_Dt   ! IAU weight for the current time step
755            !
756            IF(lwp) THEN
757               WRITE(numout,*) 
758               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
759                  &  kt,' with IAU weight = ', wgtiau(it)
760               WRITE(numout,*) '~~~~~~~~~~~~'
761            ENDIF
762            !
763            ! Save the tendency associated with the IAU weighted SSH increment
764            ! (applied in dynspg.*)
765#if defined key_asminc
766            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
767#endif
768            !
769         ELSE IF( kt == nitiaufin_r+1 ) THEN
770            !
771            ! test on ssh_bkginc needed as ssh_asm_inc is called twice by time step
772            IF (ALLOCATED(ssh_bkginc)) DEALLOCATE( ssh_bkginc )
773            !
774#if defined key_asminc
775            ssh_iau(:,:) = 0._wp
776#endif
777            !
778         ENDIF
779         !                          !-----------------------------------------
780      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
781         !                          !-----------------------------------------
782         !
783         IF ( kt == nitdin_r ) THEN
784            !
785            l_1st_euler = .TRUE.                            ! Force Euler forward step
786            !
787            ssh(:,:,Kmm) = ssh_bkg(:,:) + ssh_bkginc(:,:)   ! Initialize the now fields the background + increment
788            !
789            ssh(:,:,Kbb) = ssh(:,:,Kmm)                        ! Update before fields
790#if ! defined key_qco
791            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
792#endif
793!!gm why not e3u(:,:,:,Kbb), e3v(:,:,:,Kbb), gdept(:,:,:,Kbb) ????
794            !
795            DEALLOCATE( ssh_bkg    )
796            DEALLOCATE( ssh_bkginc )
797            !
798         ENDIF
799         !
800      ENDIF
801      !
802   END SUBROUTINE ssh_asm_inc
803
804
805   SUBROUTINE ssh_asm_div( kt, Kbb, Kmm, phdivn )
806      !!----------------------------------------------------------------------
807      !!                  ***  ROUTINE ssh_asm_div  ***
808      !!
809      !! ** Purpose :   ssh increment with z* is incorporated via a correction of the local divergence         
810      !!                across all the water column
811      !!
812      !! ** Method  :
813      !!                CAUTION : sshiau is positive (inflow) decreasing the
814      !!                          divergence and expressed in m/s
815      !!
816      !! ** Action  :   phdivn   decreased by the ssh increment
817      !!----------------------------------------------------------------------
818      INTEGER, INTENT(IN) :: kt                               ! ocean time-step index
819      INTEGER, INTENT(IN) :: Kbb, Kmm                         ! time level indices
820      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   phdivn   ! horizontal divergence
821      !!
822      INTEGER  ::   jk                                        ! dummy loop index
823      REAL(wp), DIMENSION(:,:)  , POINTER       ::   ztim     ! local array
824      !!----------------------------------------------------------------------
825      !
826#if defined key_asminc
827      CALL ssh_asm_inc( kt, Kbb, Kmm ) !==   (calculate increments)
828      !
829      IF( ln_linssh ) THEN
830         phdivn(:,:,1) = phdivn(:,:,1) - ssh_iau(:,:) / e3t(:,:,1,Kmm) * tmask(:,:,1)
831      ELSE
832         ALLOCATE( ztim(jpi,jpj) )
833         ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) )
834         DO jk = 1, jpkm1                                 
835            phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) 
836         END DO
837         !
838         DEALLOCATE(ztim)
839      ENDIF
840#endif
841      !
842   END SUBROUTINE ssh_asm_div
843
844
845   SUBROUTINE seaice_asm_inc( kt, kindic )
846      !!----------------------------------------------------------------------
847      !!                    ***  ROUTINE seaice_asm_inc  ***
848      !!         
849      !! ** Purpose : Apply the sea ice assimilation increment.
850      !!
851      !! ** Method  : Direct initialization or Incremental Analysis Updating.
852      !!
853      !! ** Action  :
854      !!
855      !!----------------------------------------------------------------------
856      INTEGER, INTENT(in)           ::   kt       ! Current time step
857      INTEGER, INTENT(in), OPTIONAL ::   kindic   ! flag for disabling the deallocation
858      !
859      INTEGER  ::   ji, jj
860      INTEGER  ::   it
861      REAL(wp) ::   zincwgt   ! IAU weight for current time step
862#if defined key_si3
863      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zofrld, zohicif, zseaicendg, zhicifinc
864      REAL(wp) ::   zhicifmin = 0.5_wp      ! ice minimum depth in metres
865#endif
866      !!----------------------------------------------------------------------
867      !
868      !                             !-----------------------------------------
869      IF ( ln_asmiau ) THEN         ! Incremental Analysis Updating
870         !                          !-----------------------------------------
871         !
872         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
873            !
874            it = kt - nit000 + 1
875            zincwgt = wgtiau(it)      ! IAU weight for the current time step
876            ! note this is not a tendency so should not be divided by rn_Dt (as with the tracer and other increments)
877            !
878            IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
879               IF(lwp) THEN
880                  WRITE(numout,*)
881                  WRITE(numout,*) 'seaice_asm_inc : sea ice conc IAU at time step = ', kt,' with IAU weight = ', wgtiau(it)
882                  WRITE(numout,*) '~~~~~~~~~~~~'
883               ENDIF
884            ENDIF
885            !
886            ! Sea-ice : SI3 case
887            !
888#if defined key_si3
889            DO_2D( 0, 0, 0, 0 )
890               zofrld (ji,jj) = 1._wp - at_i(ji,jj)
891               zohicif(ji,jj) = hm_i(ji,jj)
892               !
893               at_i  (ji,jj) = 1. - MIN( MAX( 1.-at_i  (ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp)
894               at_i_b(ji,jj) = 1. - MIN( MAX( 1.-at_i_b(ji,jj) - seaice_bkginc(ji,jj) * zincwgt, 0.0_wp), 1.0_wp)
895               fr_i(ji,jj) = at_i(ji,jj)        ! adjust ice fraction
896               !
897               zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj))   ! find out actual sea ice nudge applied
898            END_2D
899            !
900            ! Nudge sea ice depth to bring it up to a required minimum depth
901            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(A2D(0)) < zhicifmin )
902               zhicifinc(:,:) = (zhicifmin - hm_i(A2D(0))) * zincwgt
903            ELSEWHERE
904               zhicifinc(:,:) = 0.0_wp
905            END WHERE
906            !
907            ! nudge ice depth
908            DO_2D( 0, 0, 0, 0 )
909               hm_i (ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj)
910            END_2D
911            !
912            ! seaice salinity balancing (to add)
913#endif
914            !
915#if defined key_cice && defined key_asminc
916            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
917            DO_2D( 0, 0, 0, 0 )
918               ndaice_da(ji,jj) = seaice_bkginc(ji,jj) * zincwgt / rn_Dt
919            END_2D
920#endif
921            !
922            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
923               IF ( kt == nitiaufin_r ) THEN
924                  DEALLOCATE( seaice_bkginc )
925               ENDIF
926            ENDIF
927            !
928         ELSE
929            !
930#if defined key_cice && defined key_asminc
931            DO_2D( 0, 0, 0, 0 )
932               ndaice_da(ji,jj) = 0._wp        ! Sea-ice : CICE case. Zero ice increment tendency into CICE
933            END_2D
934#endif
935            !
936         ENDIF
937         !                          !-----------------------------------------
938      ELSEIF ( ln_asmdin ) THEN     ! Direct Initialization
939         !                          !-----------------------------------------
940         !
941         IF ( kt == nitdin_r ) THEN
942            !
943            l_1st_euler = .TRUE.              ! Force Euler forward step
944            !
945            ! Sea-ice : SI3 case
946            !
947#if defined key_si3
948            DO_2D( 0, 0, 0, 0 )
949               zofrld (ji,jj) = 1._wp - at_i(ji,jj)
950               zohicif(ji,jj) = hm_i(ji,jj)
951               !
952               ! Initialize the now fields the background + increment
953               at_i(ji,jj) = 1. - MIN( MAX( 1.-at_i(ji,jj) - seaice_bkginc(ji,jj), 0.0_wp), 1.0_wp)
954               at_i_b(ji,jj) = at_i(ji,jj)
955               fr_i(ji,jj) = at_i(ji,jj)        ! adjust ice fraction
956               !
957               zseaicendg(ji,jj) = zofrld(ji,jj) - (1. - at_i(ji,jj))   ! find out actual sea ice nudge applied
958            END_2D
959            !
960            ! Nudge sea ice depth to bring it up to a required minimum depth
961            WHERE( zseaicendg(:,:) > 0.0_wp .AND. hm_i(A2D(0)) < zhicifmin )
962               zhicifinc(:,:) = zhicifmin - hm_i(A2D(0))
963            ELSEWHERE
964               zhicifinc(:,:) = 0.0_wp
965            END WHERE
966            !
967            ! nudge ice depth
968            DO_2D( 0, 0, 0, 0 )
969               hm_i(ji,jj) = hm_i (ji,jj) + zhicifinc(ji,jj)
970            END_2D
971            !
972            ! seaice salinity balancing (to add)
973#endif
974            !
975#if defined key_cice && defined key_asminc
976            ! Sea-ice : CICE case. Pass ice increment tendency into CICE
977            DO_2D( 0, 0, 0, 0 )
978               ndaice_da(ji,jj) = seaice_bkginc(ji,jj) / rn_Dt
979            END_2D
980#endif
981            IF( ntile == 0 .OR. ntile == nijtile )  THEN                ! Do only on the last tile
982               IF ( .NOT. PRESENT(kindic) ) THEN
983                  DEALLOCATE( seaice_bkginc )
984               END IF
985            ENDIF
986            !
987         ELSE
988            !
989#if defined key_cice && defined key_asminc
990            DO_2D( 0, 0, 0, 0 )
991               ndaice_da(ji,jj) = 0._wp     ! Sea-ice : CICE case. Zero ice increment tendency into CICE
992            END_2D
993#endif
994            !
995         ENDIF
996
997!#if defined defined key_si3 || defined key_cice
998!
999!            IF (ln_seaicebal ) THEN       
1000!             !! balancing salinity increments
1001!             !! simple case from limflx.F90 (doesn't include a mass flux)
1002!             !! assumption is that as ice concentration is reduced or increased
1003!             !! the snow and ice depths remain constant
1004!             !! note that snow is being created where ice concentration is being increased
1005!             !! - could be more sophisticated and
1006!             !! not do this (but would need to alter h_snow)
1007!
1008!             usave(:,:,:)=sb(:,:,:)   ! use array as a temporary store
1009!
1010!             DO jj = 1, jpj
1011!               DO ji = 1, jpi
1012!           ! calculate change in ice and snow mass per unit area
1013!           ! positive values imply adding salt to the ocean (results from ice formation)
1014!           ! fwf : ice formation and melting
1015!
1016!                 zfons = ( -nfresh_da(ji,jj)*soce + nfsalt_da(ji,jj) )*rn_Dt
1017!
1018!           ! change salinity down to mixed layer depth
1019!                 mld=hmld_kara(ji,jj)
1020!
1021!           ! prevent small mld
1022!           ! less than 10m can cause salinity instability
1023!                 IF (mld < 10) mld=10
1024!
1025!           ! set to bottom of a level
1026!                 DO jk = jpk-1, 2, -1
1027!                   IF ((mld > gdepw(ji,jj,jk,Kmm)) .and. (mld < gdepw(ji,jj,jk+1,Kmm))) THEN
1028!                     mld=gdepw(ji,jj,jk+1,Kmm)
1029!                     jkmax=jk
1030!                   ENDIF
1031!                 ENDDO
1032!
1033!            ! avoid applying salinity balancing in shallow water or on land
1034!            !
1035!
1036!            ! dsal_ocn (psu kg m^-2) / (kg m^-3 * m)
1037!
1038!                 dsal_ocn=0.0_wp
1039!                 sal_thresh=5.0_wp        ! minimum salinity threshold for salinity balancing
1040!
1041!                 if (tmask(ji,jj,1) > 0 .AND. tmask(ji,jj,jkmax) > 0 ) &
1042!                              dsal_ocn = zfons / (rhop(ji,jj,1) * mld)
1043!
1044!           ! put increments in for levels in the mixed layer
1045!           ! but prevent salinity below a threshold value
1046!
1047!                   DO jk = 1, jkmax             
1048!
1049!                     IF (dsal_ocn > 0.0_wp .or. sb(ji,jj,jk)+dsal_ocn > sal_thresh) THEN
1050!                           sb(ji,jj,jk) = sb(ji,jj,jk) + dsal_ocn
1051!                           sn(ji,jj,jk) = sn(ji,jj,jk) + dsal_ocn
1052!                     ENDIF
1053!
1054!                   ENDDO
1055!
1056!      !            !  salt exchanges at the ice/ocean interface
1057!      !            zpmess         = zfons / rDt_ice    ! rDt_ice is ice timestep
1058!      !
1059!      !! Adjust fsalt. A +ve fsalt means adding salt to ocean
1060!      !!           fsalt(ji,jj) =  fsalt(ji,jj) + zpmess     ! adjust fsalt 
1061!      !!               
1062!      !!           emps(ji,jj) = emps(ji,jj) + zpmess        ! or adjust emps (see icestp1d)
1063!      !!                                                     ! E-P (kg m-2 s-2)
1064!      !            emp(ji,jj) = emp(ji,jj) + zpmess          ! E-P (kg m-2 s-2)
1065!               ENDDO !ji
1066!             ENDDO !jj!
1067!
1068!            ENDIF !ln_seaicebal
1069!
1070!#endif
1071         !
1072      ENDIF
1073      !
1074   END SUBROUTINE seaice_asm_inc
1075   
1076   !!======================================================================
1077END MODULE asminc
Note: See TracBrowser for help on using the repository browser.