source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 8029

Last change on this file since 8029 was 8029, checked in by lovato, 4 years ago

3.4 stable: apply macro selection as in 3.6_stable (#1528)

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