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

Last change on this file since 2287 was 2287, checked in by smasson, 10 years ago

update licence of all NEMO files…

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