source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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