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 trunk/NEMOGCM/NEMO/OPA_SRC/ASM – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 3604

Last change on this file since 3604 was 3604, checked in by rblod, 10 years ago

Adding routines and modules for TAM - Ticket #1005

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