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 branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 6 years ago

Remove svn keywords

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