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 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

Merge of 3.4beta into the trunk
  • Property svn:keywords set to Id
File size: 37.5 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   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   'key_asminc' : Switch on the assimilation increment interface
16   !!----------------------------------------------------------------------
17   !!   asm_inc_init : Initialize the increment arrays and IAU weights
18   !!   calc_date    : Compute the calendar date YYYYMMDD on a given step
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   !!----------------------------------------------------------------------
23   USE wrk_nemo         ! Memory Allocation
24   USE par_oce          ! Ocean space and time domain variables
25   USE dom_oce          ! Ocean space and time domain
26   USE oce              ! Dynamics and active tracers defined in memory
27   USE divcur           ! Horizontal divergence and relative vorticity
28   USE ldfdyn_oce       ! ocean dynamics: lateral physics
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
37   IMPLICIT NONE
38   PRIVATE
39   
40   PUBLIC   asm_inc_init   !: Initialize the increment arrays and IAU weights
41   PUBLIC   calc_date      !: Compute the calendar date YYYYMMDD on a given step
42   PUBLIC   tra_asm_inc    !: Apply the tracer (T and S) increments
43   PUBLIC   dyn_asm_inc    !: Apply the dynamic (u and v) increments
44   PUBLIC   ssh_asm_inc    !: Apply the SSH increment
45
46#if defined key_asminc
47    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .TRUE.   !: Logical switch for assimilation increment interface
48#else
49    LOGICAL, PUBLIC, PARAMETER :: lk_asminc = .FALSE.  !: No assimilation increments
50#endif
51   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE. !: No output of the background state fields
52   LOGICAL, PUBLIC :: ln_trjwri = .FALSE. !: No output of the state trajectory fields
53   LOGICAL, PUBLIC :: ln_asmiau = .FALSE. !: No applying forcing with an assimilation increment
54   LOGICAL, PUBLIC :: ln_asmdin = .FALSE. !: No direct initialization
55   LOGICAL, PUBLIC :: ln_trainc = .FALSE. !: No tracer (T and S) assimilation increments
56   LOGICAL, PUBLIC :: ln_dyninc = .FALSE. !: No dynamics (u and v) assimilation increments
57   LOGICAL, PUBLIC :: ln_sshinc = .FALSE. !: No sea surface height assimilation increment
58   LOGICAL, PUBLIC :: ln_salfix = .FALSE. !: Apply minimum salinity check
59   INTEGER, PUBLIC :: nn_divdmp = 0       !: Apply divergence damping filter nn_divdmp times
60
61   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkg   , s_bkg      !: Background temperature and salinity
62   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkg   , v_bkg      !: Background u- & v- velocity components
63   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   t_bkginc, s_bkginc   !: Increment to the background T & S
64   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::   u_bkginc, v_bkginc   !: Increment to the u- & v-components
65   REAL(wp), PUBLIC, DIMENSION(:)    , ALLOCATABLE ::   wgtiau               !: IAU weights for each time step
66#if defined key_asminc
67   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE ::   ssh_iau           !: IAU-weighted sea surface height increment
68#endif
69   !                                !!! time steps relative to the cycle interval [0,nitend-nit000-1]
70   INTEGER , PUBLIC ::   nitbkg      !: Time step of the background state used in the Jb term
71   INTEGER , PUBLIC ::   nitdin      !: Time step of the background state for direct initialization
72   INTEGER , PUBLIC ::   nitiaustr   !: Time step of the start of the IAU interval
73   INTEGER , PUBLIC ::   nitiaufin   !: Time step of the end of the IAU interval
74   !
75   INTEGER , PUBLIC ::   niaufn      !: Type of IAU weighing function: = 0   Constant weighting
76   !                                 !: = 1   Linear hat-like, centred in middle of IAU interval
77   REAL(wp), PUBLIC ::   salfixmin   !: Ensure that the salinity is larger than this  value if (ln_salfix)
78
79   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment
80
81   !! * Substitutions
82#  include "domzgr_substitute.h90"
83#  include "ldfdyn_substitute.h90"
84#  include "vectopt_loop_substitute.h90"
85
86   !!----------------------------------------------------------------------
87   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
88   !! $Id$
89   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
90   !!----------------------------------------------------------------------
91CONTAINS
92
93   SUBROUTINE asm_inc_init
94      !!----------------------------------------------------------------------
95      !!                    ***  ROUTINE asm_inc_init  ***
96      !!         
97      !! ** Purpose : Initialize the assimilation increment and IAU weights.
98      !!
99      !! ** Method  : Initialize the assimilation increment and IAU weights.
100      !!
101      !! ** Action  :
102      !!----------------------------------------------------------------------
103      !!
104      !!
105      INTEGER :: ji,jj,jk
106      INTEGER :: jt
107      INTEGER :: imid
108      INTEGER :: inum
109      INTEGER :: iiauper         ! Number of time steps in the IAU period
110      INTEGER :: icycper         ! Number of time steps in the cycle
111      INTEGER :: iitend_date     ! Date YYYYMMDD of final time step
112      INTEGER :: iitbkg_date     ! Date YYYYMMDD of background time step for Jb term
113      INTEGER :: iitdin_date     ! Date YYYYMMDD of background time step for DI
114      INTEGER :: iitiaustr_date  ! Date YYYYMMDD of IAU interval start time step
115      INTEGER :: iitiaufin_date  ! Date YYYYMMDD of IAU interval final time step
116
117      REAL(wp) :: znorm        ! Normalization factor for IAU weights
118      REAL(wp) :: ztotwgt      ! Value of time-integrated IAU weights
119                               ! (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), POINTER, DIMENSION(:,:) :: hdiv
126      !!
127      NAMELIST/nam_asminc/ ln_bkgwri, ln_trjwri,                           &
128         &                 ln_trainc, ln_dyninc, ln_sshinc,                &
129         &                 ln_asmdin, ln_asmiau,                           &
130         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   &
131         &                 nittrjfrq, ln_salfix, salfixmin,                &
132         &                 nn_divdmp
133      !!----------------------------------------------------------------------
134
135      !-----------------------------------------------------------------------
136      ! Read Namelist nam_asminc : assimilation increment interface
137      !-----------------------------------------------------------------------
138
139      ! Set default values
140      ln_bkgwri = .FALSE.
141      ln_trjwri = .FALSE.
142      ln_trainc = .FALSE.
143      ln_dyninc = .FALSE.
144      ln_sshinc = .FALSE.
145      ln_asmdin = .FALSE.
146      ln_asmiau = .TRUE.
147      ln_salfix = .FALSE.
148      salfixmin = -9999
149      nitbkg    = 0
150      nitdin    = 0     
151      nitiaustr = 1
152      nitiaufin = 150      ! = 10 days with ORCA2
153      niaufn    = 0
154      nittrjfrq = 1
155
156      REWIND ( numnam )
157      READ   ( numnam, nam_asminc )
158
159      ! Control print
160      IF(lwp) THEN
161         WRITE(numout,*)
162         WRITE(numout,*) 'asm_inc_init : Assimilation increment initialization :'
163         WRITE(numout,*) '~~~~~~~~~~~~'
164         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters'
165         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri
166         WRITE(numout,*) '      Logical switch for writing out state trajectory          ln_trjwri = ', ln_trjwri
167         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc
168         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc
169         WRITE(numout,*) '      Logical switch for applying SSH increments               ln_sshinc = ', ln_sshinc
170         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin
171         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau
172         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg
173         WRITE(numout,*) '      Timestep of background for DI in [0,nitend-nit000-1]     nitdin    = ', nitdin
174         WRITE(numout,*) '      Timestep of start of IAU interval in [0,nitend-nit000-1] nitiaustr = ', nitiaustr
175         WRITE(numout,*) '      Timestep of end of IAU interval in [0,nitend-nit000-1]   nitiaufin = ', nitiaufin
176         WRITE(numout,*) '      Type of IAU weighting function                           niaufn    = ', niaufn
177         WRITE(numout,*) '      Frequency of trajectory output for 4D-VAR                nittrjfrq = ', nittrjfrq
178         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix
179         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin
180      ENDIF
181
182      nitbkg_r    = nitbkg    + nit000 - 1  ! Background time referenced to nit000
183      nitdin_r    = nitdin    + nit000 - 1  ! Background time for DI referenced to nit000
184      nitiaustr_r = nitiaustr + nit000 - 1  ! Start of IAU interval referenced to nit000
185      nitiaufin_r = nitiaufin + nit000 - 1  ! End of IAU interval referenced to nit000
186
187      iiauper = nitiaufin_r - nitiaustr_r + 1  ! IAU interval length
188      icycper = nitend      - nit000      + 1  ! Cycle interval length
189
190      ! Date of final time step
191      CALL calc_date( nit000, nitend, ndate0, iitend_date )
192
193      ! Background time for Jb referenced to ndate0
194      CALL calc_date( nit000, nitbkg_r, ndate0, iitbkg_date )
195
196      ! Background time for DI referenced to ndate0
197      CALL calc_date( nit000, nitdin_r, ndate0, iitdin_date )
198
199      ! IAU start time referenced to ndate0
200      CALL calc_date( nit000, nitiaustr_r, ndate0, iitiaustr_date )
201
202      ! IAU end time referenced to ndate0
203      CALL calc_date( nit000, nitiaufin_r, ndate0, iitiaufin_date )
204
205      IF(lwp) THEN
206         WRITE(numout,*)
207         WRITE(numout,*) '   Time steps referenced to current cycle:'
208         WRITE(numout,*) '       iitrst      = ', nit000 - 1
209         WRITE(numout,*) '       nit000      = ', nit000
210         WRITE(numout,*) '       nitend      = ', nitend
211         WRITE(numout,*) '       nitbkg_r    = ', nitbkg_r
212         WRITE(numout,*) '       nitdin_r    = ', nitdin_r
213         WRITE(numout,*) '       nitiaustr_r = ', nitiaustr_r
214         WRITE(numout,*) '       nitiaufin_r = ', nitiaufin_r
215         WRITE(numout,*) '       nittrjfrq   = ', nittrjfrq
216         WRITE(numout,*)
217         WRITE(numout,*) '   Dates referenced to current cycle:'
218         WRITE(numout,*) '       ndastp         = ', ndastp
219         WRITE(numout,*) '       ndate0         = ', ndate0
220         WRITE(numout,*) '       iitend_date    = ', iitend_date
221         WRITE(numout,*) '       iitbkg_date    = ', iitbkg_date
222         WRITE(numout,*) '       iitdin_date    = ', iitdin_date
223         WRITE(numout,*) '       iitiaustr_date = ', iitiaustr_date
224         WRITE(numout,*) '       iitiaufin_date = ', iitiaufin_date
225      ENDIF
226
227      IF ( nacc /= 0 ) &
228         & CALL ctl_stop( ' nacc /= 0 and key_asminc :',  &
229         &                ' Assimilation increments have only been implemented', &
230         &                ' for synchronous time stepping' )
231
232      IF ( ( ln_asmdin ).AND.( ln_asmiau ) )   &
233         & CALL ctl_stop( ' ln_asmdin and ln_asmiau :', &
234         &                ' Choose Direct Initialization OR Incremental Analysis Updating')
235
236      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) &
237           .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) ) &
238         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc and ln_sshinc is set to .true.', &
239         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', &
240         &                ' Inconsistent options')
241
242      IF ( ( ln_bkgwri ).AND.( ( ln_asmdin ).OR.( ln_asmiau ) ) )  &
243         & CALL ctl_stop( ' ln_bkgwri and either ln_asmdin or ln_asmiau are set to .true.:', &
244         &                ' The background state must be written before applying the increments')
245
246      IF ( ( niaufn /= 0 ).AND.( niaufn /= 1 ) ) &
247         & CALL ctl_stop( ' niaufn /= 0 or niaufn /=1 :',  &
248         &                ' Type IAU weighting function is invalid')
249
250      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ) &
251         &                     )  &
252         & CALL ctl_warn( ' ln_trainc, ln_dyninc and ln_sshinc are set to .false. :', &
253         &                ' The assimilation increments are not applied')
254
255      IF ( ( ln_asmiau ).AND.( nitiaustr == nitiaufin ) ) &
256         & CALL ctl_stop( ' nitiaustr = nitiaufin :',  &
257         &                ' IAU interval is of zero length')
258
259      IF ( ( ln_asmiau ).AND.( ( nitiaustr_r < nit000 ).OR.( nitiaufin_r > nitend ) ) ) &
260         & CALL ctl_stop( ' nitiaustr or nitiaufin :',  &
261         &                ' IAU starting or final time step is outside the cycle interval', &
262         &                 ' Valid range nit000 to nitend')
263
264      IF ( ( nitbkg_r < nit000 - 1 ).OR.( nitbkg_r > nitend ) ) &
265         & CALL ctl_stop( ' nitbkg :',  &
266         &                ' Background time step is outside the cycle interval')
267
268      IF ( ( nitdin_r < nit000 - 1 ).OR.( nitdin_r > nitend ) ) &
269         & CALL ctl_stop( ' nitdin :',  &
270         &                ' Background time step for Direct Initialization is outside', &
271         &                ' the cycle interval')
272
273      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further
274
275      !--------------------------------------------------------------------
276      ! Initialize the Incremental Analysis Updating weighting function
277      !--------------------------------------------------------------------
278
279      IF ( ln_asmiau ) THEN
280
281         ALLOCATE( wgtiau( icycper ) )
282
283         wgtiau(:) = 0.0
284
285         IF ( niaufn == 0 ) THEN
286
287            !---------------------------------------------------------
288            ! Constant IAU forcing
289            !---------------------------------------------------------
290
291            DO jt = 1, iiauper
292               wgtiau(jt+nitiaustr-1) = 1.0 / REAL( iiauper )
293            END DO
294
295         ELSEIF ( niaufn == 1 ) THEN
296
297            !---------------------------------------------------------
298            ! Linear hat-like, centred in middle of IAU interval
299            !---------------------------------------------------------
300
301            ! Compute the normalization factor
302            znorm = 0.0
303            IF ( MOD( iiauper, 2 ) == 0 ) THEN  ! Even number of time steps in IAU interval
304               imid = iiauper / 2 
305               DO jt = 1, imid
306                  znorm = znorm + REAL( jt )
307               END DO
308               znorm = 2.0 * znorm
309            ELSE                               ! Odd number of time steps in IAU interval
310               imid = ( iiauper + 1 ) / 2       
311               DO jt = 1, imid - 1
312                  znorm = znorm + REAL( jt )
313               END DO
314               znorm = 2.0 * znorm + REAL( imid )
315            ENDIF
316            znorm = 1.0 / znorm
317
318            DO jt = 1, imid - 1
319               wgtiau(jt+nitiaustr-1) = REAL( jt ) * znorm
320            END DO
321            DO jt = imid, iiauper
322               wgtiau(jt+nitiaustr-1) = REAL( iiauper - jt + 1 ) * znorm
323            END DO
324
325         ENDIF
326
327         ! Test that the integral of the weights over the weighting interval equals 1
328          IF(lwp) THEN
329             WRITE(numout,*)
330             WRITE(numout,*) 'asm_inc_init : IAU weights'
331             WRITE(numout,*) '~~~~~~~~~~~~'
332             WRITE(numout,*) '             time step         IAU  weight'
333             WRITE(numout,*) '             =========     ====================='
334             ztotwgt = 0.0
335             DO jt = 1, icycper
336                ztotwgt = ztotwgt + wgtiau(jt)
337                WRITE(numout,*) '         ', jt, '       ', wgtiau(jt) 
338             END DO   
339             WRITE(numout,*) '         ==================================='
340             WRITE(numout,*) '         Time-integrated weight = ', ztotwgt
341             WRITE(numout,*) '         ==================================='
342          ENDIF
343         
344      ENDIF
345
346      !--------------------------------------------------------------------
347      ! Allocate and initialize the increment arrays
348      !--------------------------------------------------------------------
349
350      ALLOCATE( t_bkginc(jpi,jpj,jpk) )
351      ALLOCATE( s_bkginc(jpi,jpj,jpk) )
352      ALLOCATE( u_bkginc(jpi,jpj,jpk) )
353      ALLOCATE( v_bkginc(jpi,jpj,jpk) )
354      ALLOCATE( ssh_bkginc(jpi,jpj)   )
355#if defined key_asminc
356      ALLOCATE( ssh_iau(jpi,jpj)      )
357#endif
358      t_bkginc(:,:,:) = 0.0
359      s_bkginc(:,:,:) = 0.0
360      u_bkginc(:,:,:) = 0.0
361      v_bkginc(:,:,:) = 0.0
362      ssh_bkginc(:,:) = 0.0
363#if defined key_asminc
364      ssh_iau(:,:)    = 0.0
365#endif
366      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) ) THEN
367
368         !--------------------------------------------------------------------
369         ! Read the increments from file
370         !--------------------------------------------------------------------
371
372         CALL iom_open( c_asminc, inum )
373
374         CALL iom_get( inum, 'time', zdate_inc ) 
375
376         CALL iom_get( inum, 'z_inc_dateb', z_inc_dateb )
377         CALL iom_get( inum, 'z_inc_datef', z_inc_datef )
378         z_inc_dateb = zdate_inc
379         z_inc_datef = zdate_inc
380
381         IF(lwp) THEN
382            WRITE(numout,*) 
383            WRITE(numout,*) 'asm_inc_init : Assimilation increments valid ', &
384               &            ' between dates ', NINT( z_inc_dateb ),' and ',  &
385               &            NINT( z_inc_datef )
386            WRITE(numout,*) '~~~~~~~~~~~~'
387         ENDIF
388
389         IF (     ( NINT( z_inc_dateb ) < ndastp      ) &
390            & .OR.( NINT( z_inc_datef ) > iitend_date ) ) &
391            & CALL ctl_warn( ' Validity time of assimilation increments is ', &
392            &                ' outside the assimilation interval' )
393
394         IF ( ( ln_asmdin ).AND.( NINT( zdate_inc ) /= iitdin_date ) ) &
395            & CALL ctl_warn( ' Validity time of assimilation increments does ', &
396            &                ' not agree with Direct Initialization time' )
397
398         IF ( ln_trainc ) THEN   
399            CALL iom_get( inum, jpdom_autoglo, 'bckint', t_bkginc, 1 )
400            CALL iom_get( inum, jpdom_autoglo, 'bckins', s_bkginc, 1 )
401            ! Apply the masks
402            t_bkginc(:,:,:) = t_bkginc(:,:,:) * tmask(:,:,:)
403            s_bkginc(:,:,:) = s_bkginc(:,:,:) * tmask(:,:,:)
404            ! Set missing increments to 0.0 rather than 1e+20
405            ! to allow for differences in masks
406            WHERE( ABS( t_bkginc(:,:,:) ) > 1.0e+10 ) t_bkginc(:,:,:) = 0.0
407            WHERE( ABS( s_bkginc(:,:,:) ) > 1.0e+10 ) s_bkginc(:,:,:) = 0.0
408         ENDIF
409
410         IF ( ln_dyninc ) THEN   
411            CALL iom_get( inum, jpdom_autoglo, 'bckinu', u_bkginc, 1 )             
412            CALL iom_get( inum, jpdom_autoglo, 'bckinv', v_bkginc, 1 )             
413            ! Apply the masks
414            u_bkginc(:,:,:) = u_bkginc(:,:,:) * umask(:,:,:)
415            v_bkginc(:,:,:) = v_bkginc(:,:,:) * vmask(:,:,:)
416            ! Set missing increments to 0.0 rather than 1e+20
417            ! to allow for differences in masks
418            WHERE( ABS( u_bkginc(:,:,:) ) > 1.0e+10 ) u_bkginc(:,:,:) = 0.0
419            WHERE( ABS( v_bkginc(:,:,:) ) > 1.0e+10 ) v_bkginc(:,:,:) = 0.0
420         ENDIF
421       
422         IF ( ln_sshinc ) THEN
423            CALL iom_get( inum, jpdom_autoglo, 'bckineta', ssh_bkginc, 1 )
424            ! Apply the masks
425            ssh_bkginc(:,:) = ssh_bkginc(:,:) * tmask(:,:,1)
426            ! Set missing increments to 0.0 rather than 1e+20
427            ! to allow for differences in masks
428            WHERE( ABS( ssh_bkginc(:,:) ) > 1.0e+10 ) ssh_bkginc(:,:) = 0.0
429         ENDIF
430
431         CALL iom_close( inum )
432 
433      ENDIF
434
435      !-----------------------------------------------------------------------
436      ! Apply divergence damping filter
437      !-----------------------------------------------------------------------
438
439
440      IF ( ln_dyninc .AND. nn_divdmp > 0 ) THEN
441
442      CALL wrk_alloc(jpi,jpj,hdiv) 
443
444       DO  jt = 1, nn_divdmp
445
446           DO jk = 1, jpkm1
447
448                  hdiv(:,:) = 0._wp
449
450            DO jj = 2, jpjm1
451               DO ji = fs_2, fs_jpim1   ! vector opt.
452                  hdiv(ji,jj) =   &
453                     (  e2u(ji  ,jj)*fse3u(ji  ,jj,jk) * u_bkginc(ji  ,jj,jk)       &
454                      - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * u_bkginc(ji-1,jj,jk)       &
455                      + e1v(ji,jj  )*fse3v(ji,jj  ,jk) * v_bkginc(ji,jj  ,jk)       &
456                      - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * v_bkginc(ji,jj-1,jk)  )    &
457                      / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
458               END DO
459            END DO
460
461            CALL lbc_lnk( hdiv, 'T', 1. )   ! lateral boundary cond. (no sign change)
462
463            DO jj = 2, jpjm1
464               DO ji = fs_2, fs_jpim1   ! vector opt.
465                  u_bkginc(ji,jj,jk) = u_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji+1,jj)*e2t(ji+1,jj) * hdiv(ji+1,jj)   &
466                                                                  - e1t(ji  ,jj)*e2t(ji  ,jj) * hdiv(ji  ,jj) ) &
467                                                                / e1u(ji,jj) * umask(ji,jj,jk) 
468                  v_bkginc(ji,jj,jk) = v_bkginc(ji,jj,jk) + 0.2 * ( e1t(ji,jj+1)*e2t(ji,jj+1) * hdiv(ji,jj+1)   &
469                                                                  - e1t(ji,jj  )*e2t(ji,jj  ) * hdiv(ji,jj  ) ) &
470                                                                / e2v(ji,jj) * vmask(ji,jj,jk) 
471               END DO
472            END DO
473
474           END DO
475
476       END DO
477
478       CALL wrk_dealloc(jpi,jpj,hdiv) 
479
480      ENDIF
481
482
483
484      !-----------------------------------------------------------------------
485      ! Allocate and initialize the background state arrays
486      !-----------------------------------------------------------------------
487
488      IF ( ln_asmdin ) THEN
489
490         ALLOCATE( t_bkg(jpi,jpj,jpk) )
491         ALLOCATE( s_bkg(jpi,jpj,jpk) )
492         ALLOCATE( u_bkg(jpi,jpj,jpk) )
493         ALLOCATE( v_bkg(jpi,jpj,jpk) )
494         ALLOCATE( ssh_bkg(jpi,jpj)   )
495
496         t_bkg(:,:,:) = 0.0
497         s_bkg(:,:,:) = 0.0
498         u_bkg(:,:,:) = 0.0
499         v_bkg(:,:,:) = 0.0
500         ssh_bkg(:,:) = 0.0
501
502         !--------------------------------------------------------------------
503         ! Read from file the background state at analysis time
504         !--------------------------------------------------------------------
505
506         CALL iom_open( c_asmdin, inum )
507
508         CALL iom_get( inum, 'zdate', zdate_bkg ) 
509       
510         IF(lwp) THEN
511            WRITE(numout,*) 
512            WRITE(numout,*) 'asm_inc_init : Assimilation background state valid at : ', &
513               &  NINT( zdate_bkg )
514            WRITE(numout,*) '~~~~~~~~~~~~'
515         ENDIF
516
517         IF ( NINT( zdate_bkg ) /= iitdin_date ) &
518            & CALL ctl_warn( ' Validity time of assimilation background state does', &
519            &                ' not agree with Direct Initialization time' )
520
521         IF ( ln_trainc ) THEN   
522            CALL iom_get( inum, jpdom_autoglo, 'tn', t_bkg )
523            CALL iom_get( inum, jpdom_autoglo, 'sn', s_bkg )
524            t_bkg(:,:,:) = t_bkg(:,:,:) * tmask(:,:,:)
525            s_bkg(:,:,:) = s_bkg(:,:,:) * tmask(:,:,:)
526         ENDIF
527
528         IF ( ln_dyninc ) THEN   
529            CALL iom_get( inum, jpdom_autoglo, 'un', u_bkg )
530            CALL iom_get( inum, jpdom_autoglo, 'vn', v_bkg )
531            u_bkg(:,:,:) = u_bkg(:,:,:) * umask(:,:,:)
532            v_bkg(:,:,:) = v_bkg(:,:,:) * vmask(:,:,:)
533         ENDIF
534       
535         IF ( ln_sshinc ) THEN
536            CALL iom_get( inum, jpdom_autoglo, 'sshn', ssh_bkg )
537            ssh_bkg(:,:) = ssh_bkg(:,:) * tmask(:,:,1)
538         ENDIF
539
540         CALL iom_close( inum )
541
542      ENDIF
543      !
544   END SUBROUTINE asm_inc_init
545
546
547   SUBROUTINE calc_date( kit000, kt, kdate0, kdate )
548      !!----------------------------------------------------------------------
549      !!                    ***  ROUTINE calc_date  ***
550      !!         
551      !! ** Purpose : Compute the calendar date YYYYMMDD at a given time step.
552      !!
553      !! ** Method  : Compute the calendar date YYYYMMDD at a given time step.
554      !!
555      !! ** Action  :
556      !!----------------------------------------------------------------------
557      INTEGER, INTENT(IN) :: kit000  ! Initial time step
558      INTEGER, INTENT(IN) :: kt      ! Current time step referenced to kit000
559      INTEGER, INTENT(IN) :: kdate0  ! Initial date
560      INTEGER, INTENT(OUT) :: kdate  ! Current date reference to kdate0
561      !
562      INTEGER :: iyea0    ! Initial year
563      INTEGER :: imon0    ! Initial month
564      INTEGER :: iday0    ! Initial day
565      INTEGER :: iyea     ! Current year
566      INTEGER :: imon     ! Current month
567      INTEGER :: iday     ! Current day
568      INTEGER :: idaystp  ! Number of days between initial and current date
569      INTEGER :: idaycnt  ! Day counter
570
571      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
572
573      !-----------------------------------------------------------------------
574      ! Compute the calendar date YYYYMMDD
575      !-----------------------------------------------------------------------
576
577      ! Initial date
578      iyea0 =   kdate0 / 10000
579      imon0 = ( kdate0 - ( iyea0 * 10000 ) ) / 100
580      iday0 =   kdate0 - ( iyea0 * 10000 ) - ( imon0 * 100 ) 
581
582      ! Check that kt >= kit000 - 1
583      IF ( kt < kit000 - 1 ) CALL ctl_stop( ' kt must be >= kit000 - 1')
584
585      ! If kt = kit000 - 1 then set the date to the restart date
586      IF ( kt == kit000 - 1 ) THEN
587
588         kdate = ndastp
589         RETURN
590
591      ENDIF
592
593      ! Compute the number of days from the initial date
594      idaystp = INT( REAL( kt - kit000 ) * rdt / 86400. )
595   
596      iday    = iday0
597      imon    = imon0
598      iyea    = iyea0
599      idaycnt = 0
600
601      CALL calc_month_len( iyea, imonth_len )
602
603      DO WHILE ( idaycnt < idaystp )
604         iday = iday + 1
605         IF ( iday > imonth_len(imon) )  THEN
606            iday = 1
607            imon = imon + 1
608         ENDIF
609         IF ( imon > 12 ) THEN
610            imon = 1
611            iyea = iyea + 1
612            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
613         ENDIF                 
614         idaycnt = idaycnt + 1
615      END DO
616      !
617      kdate = iyea * 10000 + imon * 100 + iday
618      !
619   END SUBROUTINE
620
621
622   SUBROUTINE calc_month_len( iyear, imonth_len )
623      !!----------------------------------------------------------------------
624      !!                    ***  ROUTINE calc_month_len  ***
625      !!         
626      !! ** Purpose : Compute the number of days in a months given a year.
627      !!
628      !! ** Method  :
629      !!----------------------------------------------------------------------
630      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
631      INTEGER :: iyear         !: year
632      !!----------------------------------------------------------------------
633      !
634      ! length of the month of the current year (from nleapy, read in namelist)
635      IF ( nleapy < 2 ) THEN
636         imonth_len(:) = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
637         IF ( nleapy == 1 ) THEN   ! we are using calendar with leap years
638            IF ( MOD(iyear, 4) == 0 .AND. ( MOD(iyear, 400) == 0 .OR. MOD(iyear, 100) /= 0 ) ) THEN
639               imonth_len(2) = 29
640            ENDIF
641         ENDIF
642      ELSE
643         imonth_len(:) = nleapy   ! all months with nleapy days per year
644      ENDIF
645      !
646   END SUBROUTINE
647
648
649   SUBROUTINE tra_asm_inc( kt )
650      !!----------------------------------------------------------------------
651      !!                    ***  ROUTINE tra_asm_inc  ***
652      !!         
653      !! ** Purpose : Apply the tracer (T and S) assimilation increments
654      !!
655      !! ** Method  : Direct initialization or Incremental Analysis Updating
656      !!
657      !! ** Action  :
658      !!----------------------------------------------------------------------
659      INTEGER, INTENT(IN) :: kt               ! Current time step
660      !
661      INTEGER :: ji,jj,jk
662      INTEGER :: it
663      REAL(wp) :: zincwgt  ! IAU weight for current time step
664      !!----------------------------------------------------------------------
665
666      IF ( ln_asmiau ) THEN
667
668         !--------------------------------------------------------------------
669         ! Incremental Analysis Updating
670         !--------------------------------------------------------------------
671
672         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
673
674            it = kt - nit000 + 1
675            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
676
677            IF(lwp) THEN
678               WRITE(numout,*) 
679               WRITE(numout,*) 'tra_asm_inc : Tracer IAU at time step = ', &
680                  &  kt,' with IAU weight = ', wgtiau(it)
681               WRITE(numout,*) '~~~~~~~~~~~~'
682            ENDIF
683
684            ! Update the tracer tendencies
685            DO jk = 1, jpkm1
686               tsa(:,:,jk,jp_tem) = tsa(:,:,jk,jp_tem) + t_bkginc(:,:,jk) * zincwgt 
687               tsa(:,:,jk,jp_sal) = tsa(:,:,jk,jp_sal) + s_bkginc(:,:,jk) * zincwgt
688            END DO
689
690            ! Salinity fix
691            IF (ln_salfix) THEN
692               DO jk = 1, jpkm1
693                  DO jj = 1, jpj
694                     DO ji= 1, jpi
695                        tsa(ji,jj,jk,jp_sal) = MAX( tsa(ji,jj,jk,jp_sal), salfixmin )
696                     END DO
697                  END DO
698               END DO
699            ENDIF
700
701         ENDIF
702
703         IF ( kt == nitiaufin_r + 1  ) THEN   ! For bias crcn to work
704            DEALLOCATE( t_bkginc )
705            DEALLOCATE( s_bkginc )
706         ENDIF
707
708
709      ELSEIF ( ln_asmdin ) THEN
710
711         !--------------------------------------------------------------------
712         ! Direct Initialization
713         !--------------------------------------------------------------------
714           
715         IF ( kt == nitdin_r ) THEN
716
717            neuler = 0  ! Force Euler forward step
718
719            ! Initialize the now fields with the background + increment
720            tsn(:,:,:,jp_tem) = t_bkg(:,:,:) + t_bkginc(:,:,:)   
721            tsn(:,:,:,jp_sal) = s_bkg(:,:,:) + s_bkginc(:,:,:)   
722
723            ! Optional salinity fix
724            IF (ln_salfix) THEN
725               DO jk = 1, jpkm1
726                  DO jj = 1, jpj
727                     DO ji= 1, jpi
728                        tsn(ji,jj,jk,jp_sal) = MAX( tsn(ji,jj,jk,jp_sal), salfixmin )
729                     END DO
730                  END DO
731               END DO
732            ENDIF
733
734            tsb(:,:,:,:) = tsn(:,:,:,:)                        ! Update before fields
735
736            CALL eos( tsb, rhd, rhop )                ! Before potential and in situ densities
737         
738            IF( ln_zps .AND. .NOT. lk_c1d ) &
739               &  CALL zps_hde( nit000, jpts, tsb,   &  ! Partial steps: before horizontal derivative
740               &                gtsu, gtsv, rhd,        &  ! of T, S, rd at the bottom ocean level
741               &                gru , grv )
742
743            DEALLOCATE( t_bkginc )
744            DEALLOCATE( s_bkginc )
745            DEALLOCATE( t_bkg    )
746            DEALLOCATE( s_bkg    )
747         ENDIF
748         ! 
749      ENDIF
750      !
751   END SUBROUTINE tra_asm_inc
752
753
754   SUBROUTINE dyn_asm_inc( kt )
755      !!----------------------------------------------------------------------
756      !!                    ***  ROUTINE dyn_asm_inc  ***
757      !!         
758      !! ** Purpose : Apply the dynamics (u and v) assimilation increments.
759      !!
760      !! ** Method  : Direct initialization or Incremental Analysis Updating.
761      !!
762      !! ** Action  :
763      !!----------------------------------------------------------------------
764      INTEGER, INTENT(IN) :: kt   ! Current time step
765      !
766      INTEGER :: jk
767      INTEGER :: it
768      REAL(wp) :: zincwgt  ! IAU weight for current time step
769      !!----------------------------------------------------------------------
770
771      IF ( ln_asmiau ) THEN
772
773         !--------------------------------------------------------------------
774         ! Incremental Analysis Updating
775         !--------------------------------------------------------------------
776
777         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
778
779            it = kt - nit000 + 1
780            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
781
782            IF(lwp) THEN
783               WRITE(numout,*) 
784               WRITE(numout,*) 'dyn_asm_inc : Dynamics IAU at time step = ', &
785                  &  kt,' with IAU weight = ', wgtiau(it)
786               WRITE(numout,*) '~~~~~~~~~~~~'
787            ENDIF
788
789            ! Update the dynamic tendencies
790            DO jk = 1, jpkm1
791               ua(:,:,jk) = ua(:,:,jk) + u_bkginc(:,:,jk) * zincwgt
792               va(:,:,jk) = va(:,:,jk) + v_bkginc(:,:,jk) * zincwgt
793            END DO
794           
795            IF ( kt == nitiaufin_r ) THEN
796               DEALLOCATE( u_bkginc )
797               DEALLOCATE( v_bkginc )
798            ENDIF
799
800         ENDIF
801
802      ELSEIF ( ln_asmdin ) THEN 
803
804         !--------------------------------------------------------------------
805         ! Direct Initialization
806         !--------------------------------------------------------------------
807         
808         IF ( kt == nitdin_r ) THEN
809
810            neuler = 0                    ! Force Euler forward step
811
812            ! Initialize the now fields with the background + increment
813            un(:,:,:) = u_bkg(:,:,:) + u_bkginc(:,:,:)
814            vn(:,:,:) = v_bkg(:,:,:) + v_bkginc(:,:,:) 
815
816            ub(:,:,:) = un(:,:,:)         ! Update before fields
817            vb(:,:,:) = vn(:,:,:)
818 
819            CALL div_cur( kt )            ! Compute divergence and curl for now fields
820
821            rotb (:,:,:) = rotn (:,:,:)   ! Update before fields
822            hdivb(:,:,:) = hdivn(:,:,:)
823
824            DEALLOCATE( u_bkg    )
825            DEALLOCATE( v_bkg    )
826            DEALLOCATE( u_bkginc )
827            DEALLOCATE( v_bkginc )
828         ENDIF
829         !
830      ENDIF
831      !
832   END SUBROUTINE dyn_asm_inc
833
834
835   SUBROUTINE ssh_asm_inc( kt )
836      !!----------------------------------------------------------------------
837      !!                    ***  ROUTINE ssh_asm_inc  ***
838      !!         
839      !! ** Purpose : Apply the sea surface height assimilation increment.
840      !!
841      !! ** Method  : Direct initialization or Incremental Analysis Updating.
842      !!
843      !! ** Action  :
844      !!----------------------------------------------------------------------
845      INTEGER, INTENT(IN) :: kt   ! Current time step
846      !
847      INTEGER :: it
848      REAL(wp) :: zincwgt  ! IAU weight for current time step
849      !!----------------------------------------------------------------------
850
851      IF ( ln_asmiau ) THEN
852
853         !--------------------------------------------------------------------
854         ! Incremental Analysis Updating
855         !--------------------------------------------------------------------
856
857         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN
858
859            it = kt - nit000 + 1
860            zincwgt = wgtiau(it) / rdt   ! IAU weight for the current time step
861
862            IF(lwp) THEN
863               WRITE(numout,*) 
864               WRITE(numout,*) 'ssh_asm_inc : SSH IAU at time step = ', &
865                  &  kt,' with IAU weight = ', wgtiau(it)
866               WRITE(numout,*) '~~~~~~~~~~~~'
867            ENDIF
868
869            ! Save the tendency associated with the IAU weighted SSH increment
870            ! (applied in dynspg.*)
871#if defined key_asminc
872            ssh_iau(:,:) = ssh_bkginc(:,:) * zincwgt
873#endif
874            IF ( kt == nitiaufin_r ) THEN
875               DEALLOCATE( ssh_bkginc )
876            ENDIF
877
878         ENDIF
879
880      ELSEIF ( ln_asmdin ) THEN
881
882         !--------------------------------------------------------------------
883         ! Direct Initialization
884         !--------------------------------------------------------------------
885
886         IF ( kt == nitdin_r ) THEN
887
888            neuler = 0                    ! Force Euler forward step
889
890            ! Initialize the now fields the background + increment
891            sshn(:,:) = ssh_bkg(:,:) + ssh_bkginc(:,:) 
892
893            sshb(:,:) = sshn(:,:)         ! Update before fields
894
895            DEALLOCATE( ssh_bkg    )
896            DEALLOCATE( ssh_bkginc )
897
898         ENDIF
899         !
900      ENDIF
901      !
902   END SUBROUTINE ssh_asm_inc
903
904   !!======================================================================
905END MODULE asminc
Note: See TracBrowser for help on using the repository browser.