source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 10 years ago

v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections

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