source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 5034

Last change on this file since 5034 was 5034, checked in by andrewryan, 6 years ago

merge with trunk

  • Property svn:keywords set to Id
File size: 60.0 KB
Line 
1MODULE diaobs
2   !!======================================================================
3   !!                       ***  MODULE diaobs  ***
4   !! Observation diagnostics: Computation of the misfit between data and
5   !!                          their model equivalent
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   'key_diaobs' : Switch on the observation diagnostic computation
10   !!----------------------------------------------------------------------
11   !!   dia_obs_init : Reading and prepare observations
12   !!   dia_obs      : Compute model equivalent to observations
13   !!   dia_obs_wri  : Write observational diagnostics
14   !!   ini_date     : Compute the initial date YYYYMMDD.HHMMSS
15   !!   fin_date     : Compute the final date YYYYMMDD.HHMMSS
16   !!----------------------------------------------------------------------
17   !! * Modules used   
18   USE wrk_nemo                 ! Memory Allocation
19   USE par_kind                 ! Precision variables
20   USE in_out_manager           ! I/O manager
21   USE par_oce
22   USE dom_oce                  ! Ocean space and time domain variables
23   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch
24   USE obs_read_prof            ! Reading and allocation of observations (Coriolis)
25   USE obs_read_sla             ! Reading and allocation of SLA observations 
26   USE obs_read_sst             ! Reading and allocation of SST observations 
27   USE obs_readmdt              ! Reading and allocation of MDT for SLA.
28   USE obs_read_seaice          ! Reading and allocation of Sea Ice observations 
29   USE obs_read_vel             ! Reading and allocation of velocity component observations
30   USE obs_prep                 ! Preparation of obs. (grid search etc).
31   USE obs_oper                 ! Observation operators
32   USE obs_write                ! Writing of observation related diagnostics
33   USE obs_grid                 ! Grid searching
34   USE obs_read_altbias         ! Bias treatment for altimeter
35   USE obs_profiles_def         ! Profile data definitions
36   USE obs_profiles             ! Profile data storage
37   USE obs_surf_def             ! Surface data definitions
38   USE obs_sla                  ! SLA data storage
39   USE obs_sst                  ! SST data storage
40   USE obs_seaice               ! Sea Ice data storage
41   USE obs_types                ! Definitions for observation types
42   USE mpp_map                  ! MPP mapping
43   USE lib_mpp                  ! For ctl_warn/stop
44
45   IMPLICIT NONE
46
47   !! * Routine accessibility
48   PRIVATE
49   PUBLIC dia_obs_init, &  ! Initialize and read observations
50      &   dia_obs,      &  ! Compute model equivalent to observations
51      &   dia_obs_wri,  &  ! Write model equivalent to observations
52      &   dia_obs_dealloc  ! Deallocate dia_obs data
53
54   !! * Shared Module variables
55   LOGICAL, PUBLIC, PARAMETER :: &
56#if defined key_diaobs
57      & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics
58#else
59      & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics
60#endif
61
62   !! * Module variables
63   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles
64   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles
65   LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set
66   LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set
67   LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles
68   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies
69   LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files
70   LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files
71   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature
72   LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature
73   LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data
74   LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files
75   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration
76   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations
77   LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data
78   LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data
79   LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data
80   LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data
81   LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files
82   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
83   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
84   LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations
85   LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land
86   LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias 
87   LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files
88   LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations
89
90   REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS
91   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS
92 
93   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method
94   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method
95
96   INTEGER, DIMENSION(imaxavtypes) :: &
97      & endailyavtypes !: ENACT data types which are daily average
98
99   INTEGER, PARAMETER :: MaxNumFiles = 1000
100   LOGICAL, DIMENSION(MaxNumFiles) :: &
101      & ln_profb_ena, & !: Is the feedback files from ENACT data ?
102   !                    !: If so use endailyavtypes
103      & ln_profb_enatim !: Change tim for 820 enact data set.
104   
105   LOGICAL, DIMENSION(MaxNumFiles) :: &
106      & ln_velfb_av   !: Is the velocity feedback files daily average?
107   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
108      & ld_enact     !: Profile data is ENACT so use endailyavtypes
109   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
110      & ld_velav     !: Velocity data is daily averaged
111   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
112      & ld_sstnight  !: SST observation corresponds to night mean
113
114   !!----------------------------------------------------------------------
115   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
116   !! $Id$
117   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
118   !!----------------------------------------------------------------------
119
120CONTAINS
121
122   SUBROUTINE dia_obs_init
123      !!----------------------------------------------------------------------
124      !!                    ***  ROUTINE dia_obs_init  ***
125      !!         
126      !! ** Purpose : Initialize and read observations
127      !!
128      !! ** Method  : Read the namelist and call reading routines
129      !!
130      !! ** Action  : Read the namelist and call reading routines
131      !!
132      !! History :
133      !!        !  06-03  (K. Mogensen) Original code
134      !!        !  06-05  (A. Weaver) Reformatted
135      !!        !  06-10  (A. Weaver) Cleaning and add controls
136      !!        !  07-03  (K. Mogensen) General handling of profiles
137      !!----------------------------------------------------------------------
138
139      IMPLICIT NONE
140
141      !! * Local declarations
142      CHARACTER(len=128) :: enactfiles(MaxNumFiles)
143      CHARACTER(len=128) :: coriofiles(MaxNumFiles)
144      CHARACTER(len=128) :: profbfiles(MaxNumFiles)
145      CHARACTER(len=128) :: sstfiles(MaxNumFiles)     
146      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles) 
147      CHARACTER(len=128) :: slafilesact(MaxNumFiles)     
148      CHARACTER(len=128) :: slafilespas(MaxNumFiles)     
149      CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
150      CHARACTER(len=128) :: seaicefiles(MaxNumFiles)           
151      CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 
152      CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)   
153      CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
154      CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
155      CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
156      CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
157      CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
158      CHARACTER(LEN=128) :: reysstname
159      CHARACTER(LEN=12)  :: reysstfmt
160      CHARACTER(LEN=128) :: bias_file
161      CHARACTER(LEN=20)  :: datestr=" ", timestr=" "
162      NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       &
163         &            ln_sla, ln_sladt, ln_slafb,                     &
164         &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       &
165         &            enactfiles, coriofiles, profbfiles,             &
166         &            slafilesact, slafilespas, slafbfiles,           &
167         &            sstfiles, sstfbfiles,                           &
168         &            ln_seaice, seaicefiles,                         &
169         &            dobsini, dobsend, n1dint, n2dint,               &
170         &            nmsshc, mdtcorr, mdtcutoff,                     &
171         &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    &
172         &            ln_sstnight,                                    &
173         &            ln_grid_search_lookup,                          &
174         &            grid_search_file, grid_search_res,              &
175         &            ln_grid_global, bias_file, ln_altbias,          &
176         &            endailyavtypes, ln_s_at_t, ln_profb_ena,        &
177         &            ln_vel3d, ln_velavcur, velavcurfiles,           &
178         &            ln_velhrcur, velhrcurfiles,                     &
179         &            ln_velavadcp, velavadcpfiles,                   &
180         &            ln_velhradcp, velhradcpfiles,                   &
181         &            ln_velfb, velfbfiles, ln_velfb_av,              &
182         &            ln_profb_enatim, ln_ignmis, ln_cl4
183
184      INTEGER :: jprofset
185      INTEGER :: jveloset
186      INTEGER :: jvar
187      INTEGER :: jnumenact
188      INTEGER :: jnumcorio
189      INTEGER :: jnumprofb
190      INTEGER :: jnumslaact
191      INTEGER :: jnumslapas
192      INTEGER :: jnumslafb
193      INTEGER :: jnumsst
194      INTEGER :: jnumsstfb
195      INTEGER :: jnumseaice
196      INTEGER :: jnumvelavcur
197      INTEGER :: jnumvelhrcur 
198      INTEGER :: jnumvelavadcp
199      INTEGER :: jnumvelhradcp   
200      INTEGER :: jnumvelfb
201      INTEGER :: ji
202      INTEGER :: jset
203      INTEGER :: ios                 ! Local integer output status for namelist read
204      LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
205
206      !-----------------------------------------------------------------------
207      ! Read namelist parameters
208      !-----------------------------------------------------------------------
209
210      enactfiles(:) = ''
211      coriofiles(:) = ''
212      profbfiles(:) = ''
213      slafilesact(:) = ''
214      slafilespas(:) = ''
215      slafbfiles(:) = ''
216      sstfiles(:)   = ''
217      sstfbfiles(:) = ''
218      seaicefiles(:) = ''
219      velcurfiles(:) = ''
220      veladcpfiles(:) = ''
221      velavcurfiles(:) = ''
222      velhrcurfiles(:) = ''
223      velavadcpfiles(:) = ''
224      velhradcpfiles(:) = ''
225      velfbfiles(:) = ''
226      velcurfiles(:) = ''
227      veladcpfiles(:) = ''
228      endailyavtypes(:) = -1
229      endailyavtypes(1) = 820
230      ln_profb_ena(:) = .FALSE.
231      ln_profb_enatim(:) = .TRUE.
232      ln_velfb_av(:) = .FALSE.
233      ln_ignmis = .FALSE.
234     
235      CALL ini_date( dobsini )
236      CALL fin_date( dobsend )
237 
238      ! Read Namelist namobs : control observation diagnostics
239      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation
240      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
241901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
242
243      REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation
244      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
245902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
246      IF(lwm) WRITE ( numond, namobs )
247
248      ! Count number of files for each type
249      IF (ln_ena) THEN
250         lmask(:) = .FALSE.
251         WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
252         jnumenact = COUNT(lmask)
253      ENDIF
254      IF (ln_cor) THEN
255         lmask(:) = .FALSE.
256         WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
257         jnumcorio = COUNT(lmask)
258      ENDIF
259      IF (ln_profb) THEN
260         lmask(:) = .FALSE.
261         WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
262         jnumprofb = COUNT(lmask)
263      ENDIF
264      IF (ln_sladt) THEN
265         lmask(:) = .FALSE.
266         WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
267         jnumslaact = COUNT(lmask)
268         lmask(:) = .FALSE.
269         WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
270         jnumslapas = COUNT(lmask)
271      ENDIF
272      IF (ln_slafb) THEN
273         lmask(:) = .FALSE.
274         WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
275         jnumslafb = COUNT(lmask)
276         lmask(:) = .FALSE.
277      ENDIF
278      IF (ln_ghrsst) THEN
279         lmask(:) = .FALSE.
280         WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
281         jnumsst = COUNT(lmask)
282      ENDIF     
283      IF (ln_sstfb) THEN
284         lmask(:) = .FALSE.
285         WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
286         jnumsstfb = COUNT(lmask)
287         lmask(:) = .FALSE.
288      ENDIF
289      IF (ln_seaice) THEN
290         lmask(:) = .FALSE.
291         WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
292         jnumseaice = COUNT(lmask)
293      ENDIF
294      IF (ln_velavcur) THEN
295         lmask(:) = .FALSE.
296         WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
297         jnumvelavcur = COUNT(lmask)
298      ENDIF
299      IF (ln_velhrcur) THEN
300         lmask(:) = .FALSE.
301         WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
302         jnumvelhrcur = COUNT(lmask)
303      ENDIF
304      IF (ln_velavadcp) THEN
305         lmask(:) = .FALSE.
306         WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
307         jnumvelavadcp = COUNT(lmask)
308      ENDIF
309      IF (ln_velhradcp) THEN
310         lmask(:) = .FALSE.
311         WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
312         jnumvelhradcp = COUNT(lmask)
313      ENDIF
314      IF (ln_velfb) THEN
315         lmask(:) = .FALSE.
316         WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
317         jnumvelfb = COUNT(lmask)
318         lmask(:) = .FALSE.
319      ENDIF
320     
321      ! Control print
322      IF(lwp) THEN
323         WRITE(numout,*)
324         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
325         WRITE(numout,*) '~~~~~~~~~~~~'
326         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
327         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d
328         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d
329         WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena
330         WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor
331         WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb
332         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla
333         WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt
334         WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb
335         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh
336         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst
337         WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst   
338         WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst
339         WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb
340         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight
341         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss
342         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice
343         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d
344         WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur
345         WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur
346         WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp
347         WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp
348         WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb
349         WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
350         WRITE(numout,*) &
351   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup
352         IF (ln_grid_search_lookup) &
353            WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file
354         IF (ln_ena) THEN
355            DO ji = 1, jnumenact
356               WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', &
357                  TRIM(enactfiles(ji))
358            END DO
359         ENDIF
360         IF (ln_cor) THEN
361            DO ji = 1, jnumcorio
362               WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', &
363                  TRIM(coriofiles(ji))
364            END DO
365         ENDIF
366         IF (ln_profb) THEN
367            DO ji = 1, jnumprofb
368               IF (ln_profb_ena(ji)) THEN
369                  WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', &
370                     TRIM(profbfiles(ji))
371               ELSE
372                  WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', &
373                     TRIM(profbfiles(ji))
374               ENDIF
375               WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji)
376            END DO
377         ENDIF
378         IF (ln_sladt) THEN
379            DO ji = 1, jnumslaact
380               WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', &
381                  TRIM(slafilesact(ji))
382            END DO
383            DO ji = 1, jnumslapas
384               WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', &
385                  TRIM(slafilespas(ji))
386            END DO
387         ENDIF
388         IF (ln_slafb) THEN
389            DO ji = 1, jnumslafb
390               WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', &
391                  TRIM(slafbfiles(ji))
392            END DO
393         ENDIF
394         IF (ln_ghrsst) THEN
395            DO ji = 1, jnumsst
396               WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', &
397                  TRIM(sstfiles(ji))
398            END DO
399         ENDIF
400         IF (ln_sstfb) THEN
401            DO ji = 1, jnumsstfb
402               WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', &
403                  TRIM(sstfbfiles(ji))
404            END DO
405         ENDIF
406         IF (ln_seaice) THEN
407            DO ji = 1, jnumseaice
408               WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', &
409                  TRIM(seaicefiles(ji))
410            END DO
411         ENDIF
412         IF (ln_velavcur) THEN
413            DO ji = 1, jnumvelavcur
414               WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', &
415                  TRIM(velavcurfiles(ji))
416            END DO
417         ENDIF
418         IF (ln_velhrcur) THEN
419            DO ji = 1, jnumvelhrcur
420               WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', &
421                  TRIM(velhrcurfiles(ji))
422            END DO
423         ENDIF
424         IF (ln_velavadcp) THEN
425            DO ji = 1, jnumvelavadcp
426               WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', &
427                  TRIM(velavadcpfiles(ji))
428            END DO
429         ENDIF
430         IF (ln_velhradcp) THEN
431            DO ji = 1, jnumvelhradcp
432               WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', &
433                  TRIM(velhradcpfiles(ji))
434            END DO
435         ENDIF
436         IF (ln_velfb) THEN
437            DO ji = 1, jnumvelfb
438               IF (ln_velfb_av(ji)) THEN
439                  WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', &
440                     TRIM(velfbfiles(ji))
441               ELSE
442                  WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', &
443                     TRIM(velfbfiles(ji))
444               ENDIF
445            END DO
446         ENDIF
447         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini
448         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend
449         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint
450         WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint
451         WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea
452         WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc
453         WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr
454         WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff
455         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias
456         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis
457         WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes
458
459      ENDIF
460     
461      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
462         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
463         RETURN
464      ENDIF
465
466      CALL obs_typ_init
467     
468      CALL mppmap_init
469     
470      ! Parameter control
471#if defined key_diaobs
472      IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
473         & ( .NOT. ln_vel3d ).AND.                                         &
474         & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
475         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ) ) THEN
476         IF(lwp) WRITE(numout,cform_war)
477         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
478            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d are all set to .FALSE.'
479         nwarn = nwarn + 1
480      ENDIF
481#endif
482
483      CALL obs_grid_setup( )
484      IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
485         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
486            &                    ' is not available')
487      ENDIF
488      IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
489         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
490            &                    ' is not available')
491      ENDIF
492
493      !-----------------------------------------------------------------------
494      ! Depending on switches read the various observation types
495      !-----------------------------------------------------------------------
496      !  - Temperature/salinity profiles
497
498      IF ( ln_t3d .OR. ln_s3d ) THEN
499
500         ! Set the number of variables for profiles to 2 (T and S)
501         nprofvars = 2
502         ! Set the number of extra variables for profiles to 1 (insitu temp).
503         nprofextr = 1
504
505         ! Count how may insitu data sets we have and allocate data.
506         jprofset = 0
507         IF ( ln_ena ) jprofset = jprofset + 1
508         IF ( ln_cor ) jprofset = jprofset + 1
509         IF ( ln_profb ) jprofset = jprofset + jnumprofb
510         nprofsets = jprofset
511         IF ( nprofsets > 0 ) THEN
512            ALLOCATE(ld_enact(nprofsets))
513            ALLOCATE(profdata(nprofsets))
514            ALLOCATE(prodatqc(nprofsets))
515         ENDIF
516
517         jprofset = 0
518         
519         ! ENACT insitu data
520
521         IF ( ln_ena ) THEN
522
523            jprofset = jprofset + 1
524           
525            ld_enact(jprofset) = .TRUE.
526
527            CALL obs_rea_pro_dri( 1, profdata(jprofset),          &
528               &                  jnumenact, enactfiles(1:jnumenact), &
529               &                  nprofvars, nprofextr,        &
530               &                  nitend-nit000+2,             &
531               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
532               &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
533               &                  kdailyavtypes = endailyavtypes )
534
535            DO jvar = 1, 2
536
537               CALL obs_prof_staend( profdata(jprofset), jvar )
538
539            END DO
540
541            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
542               &              ln_t3d, ln_s3d, ln_nea, &
543               &              kdailyavtypes=endailyavtypes )
544           
545         ENDIF
546
547         ! Coriolis insitu data
548
549         IF ( ln_cor ) THEN
550           
551            jprofset = jprofset + 1
552
553            ld_enact(jprofset) = .FALSE.
554
555            CALL obs_rea_pro_dri( 2, profdata(jprofset),          &
556               &                  jnumcorio, coriofiles(1:jnumcorio), &
557               &                  nprofvars, nprofextr,        &
558               &                  nitend-nit000+2,             &
559               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
560               &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
561
562            DO jvar = 1, 2
563
564               CALL obs_prof_staend( profdata(jprofset), jvar )
565
566            END DO
567
568            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
569                 &            ln_t3d, ln_s3d, ln_nea )
570           
571         ENDIF
572 
573         ! Feedback insitu data
574
575         IF ( ln_profb ) THEN
576           
577            DO jset = 1, jnumprofb
578               
579               jprofset = jprofset + 1
580               ld_enact (jprofset) = ln_profb_ena(jset)
581
582               CALL obs_rea_pro_dri( 0, profdata(jprofset),          &
583                  &                  1, profbfiles(jset:jset), &
584                  &                  nprofvars, nprofextr,        &
585                  &                  nitend-nit000+2,             &
586                  &                  dobsini, dobsend, ln_t3d, ln_s3d, &
587                  &                  ln_ignmis, ln_s_at_t, &
588                  &                  ld_enact(jprofset).AND.&
589                  &                  ln_profb_enatim(jset), &
590                  &                  .FALSE., kdailyavtypes = endailyavtypes )
591               
592               DO jvar = 1, 2
593                 
594                  CALL obs_prof_staend( profdata(jprofset), jvar )
595                 
596               END DO
597               
598               IF ( ld_enact(jprofset) ) THEN
599                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
600                     &              ln_t3d, ln_s3d, ln_nea, &
601                     &              kdailyavtypes = endailyavtypes )
602               ELSE
603                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
604                     &              ln_t3d, ln_s3d, ln_nea )
605               ENDIF
606               
607            END DO
608
609         ENDIF
610
611      ENDIF
612
613      !  - Sea level anomalies
614      IF ( ln_sla ) THEN
615        ! Set the number of variables for sla to 1
616         nslavars = 1
617
618         ! Set the number of extra variables for sla to 2
619         nslaextr = 2
620         
621         ! Set the number of sla data sets to 2
622         nslasets = 0
623         IF ( ln_sladt ) THEN
624            nslasets = nslasets + 2
625         ENDIF
626         IF ( ln_slafb ) THEN
627            nslasets = nslasets + jnumslafb
628         ENDIF
629         
630         ALLOCATE(sladata(nslasets))
631         ALLOCATE(sladatqc(nslasets))
632         sladata(:)%nsurf=0
633         sladatqc(:)%nsurf=0
634
635         nslasets = 0
636
637         ! AVISO SLA data
638
639         IF ( ln_sladt ) THEN
640
641            ! Active SLA observations
642           
643            nslasets = nslasets + 1
644           
645            CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
646               &              slafilesact(1:jnumslaact), &
647               &              nslavars, nslaextr, nitend-nit000+2, &
648               &              dobsini, dobsend, ln_ignmis, .FALSE. )
649            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
650               &              ln_sla, ln_nea )
651           
652            ! Passive SLA observations
653           
654            nslasets = nslasets + 1
655           
656            CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
657               &              slafilespas(1:jnumslapas), &
658               &              nslavars, nslaextr, nitend-nit000+2, &
659               &              dobsini, dobsend, ln_ignmis, .FALSE. )
660           
661            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
662               &              ln_sla, ln_nea )
663
664         ENDIF
665         
666         ! Feedback SLA data
667
668         IF ( ln_slafb ) THEN
669
670            DO jset = 1, jnumslafb
671           
672               nslasets = nslasets + 1
673           
674               CALL obs_rea_sla( 0, sladata(nslasets), 1, &
675                  &              slafbfiles(jset:jset), &
676                  &              nslavars, nslaextr, nitend-nit000+2, &
677                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
678               CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
679                  &              ln_sla, ln_nea )
680
681            END DO               
682
683         ENDIF
684         
685         CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
686           
687         ! read in altimeter bias
688         
689         IF ( ln_altbias ) THEN     
690            CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
691         ENDIF
692     
693      ENDIF
694
695      !  - Sea surface height
696      IF ( ln_ssh ) THEN
697         IF(lwp) WRITE(numout,*) ' SSH currently not available'
698      ENDIF
699
700      !  - Sea surface temperature
701      IF ( ln_sst ) THEN
702
703         ! Set the number of variables for sst to 1
704         nsstvars = 1
705
706         ! Set the number of extra variables for sst to 0
707         nsstextr = 0
708
709         nsstsets = 0
710
711         IF (ln_reysst) nsstsets = nsstsets + 1
712         IF (ln_ghrsst) nsstsets = nsstsets + 1
713         IF ( ln_sstfb ) THEN
714            nsstsets = nsstsets + jnumsstfb
715         ENDIF
716
717         ALLOCATE(sstdata(nsstsets))
718         ALLOCATE(sstdatqc(nsstsets))
719         ALLOCATE(ld_sstnight(nsstsets))
720         sstdata(:)%nsurf=0
721         sstdatqc(:)%nsurf=0   
722         ld_sstnight(:)=.false.
723
724         nsstsets = 0
725
726         IF (ln_reysst) THEN
727
728            nsstsets = nsstsets + 1
729
730            ld_sstnight(nsstsets) = ln_sstnight
731
732            CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
733               &                  nsstvars, nsstextr, &
734               &                  nitend-nit000+2, dobsini, dobsend )
735            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
736               &              ln_nea )
737
738        ENDIF
739       
740        IF (ln_ghrsst) THEN
741       
742            nsstsets = nsstsets + 1
743
744            ld_sstnight(nsstsets) = ln_sstnight
745         
746            CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
747               &              sstfiles(1:jnumsst), &
748               &              nsstvars, nsstextr, nitend-nit000+2, &
749               &              dobsini, dobsend, ln_ignmis, .FALSE. )
750            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
751               &              ln_nea )
752
753        ENDIF
754               
755         ! Feedback SST data
756
757         IF ( ln_sstfb ) THEN
758
759            DO jset = 1, jnumsstfb
760           
761               nsstsets = nsstsets + 1
762
763               ld_sstnight(nsstsets) = ln_sstnight
764           
765               CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
766                  &              sstfbfiles(jset:jset), &
767                  &              nsstvars, nsstextr, nitend-nit000+2, &
768                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
769               CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
770                  &              ln_sst, ln_nea )
771
772            END DO               
773
774         ENDIF
775
776      ENDIF
777
778      !  - Sea surface salinity
779      IF ( ln_sss ) THEN
780         IF(lwp) WRITE(numout,*) ' SSS currently not available'
781      ENDIF
782
783      !  - Sea Ice Concentration
784     
785      IF ( ln_seaice ) THEN
786
787         ! Set the number of variables for seaice to 1
788         nseaicevars = 1
789
790         ! Set the number of extra variables for seaice to 0
791         nseaiceextr = 0
792         
793         ! Set the number of data sets to 1
794         nseaicesets = 1
795
796         ALLOCATE(seaicedata(nseaicesets))
797         ALLOCATE(seaicedatqc(nseaicesets))
798         seaicedata(:)%nsurf=0
799         seaicedatqc(:)%nsurf=0
800
801         CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
802            &                 seaicefiles(1:jnumseaice), &
803            &                 nseaicevars, nseaiceextr, nitend-nit000+2, &
804            &                 dobsini, dobsend, ln_ignmis, .FALSE. )
805
806         CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
807            &                 ln_seaice, ln_nea )
808 
809      ENDIF
810
811      IF (ln_vel3d) THEN
812
813         ! Set the number of variables for profiles to 2 (U and V)
814         nvelovars = 2
815
816         ! Set the number of extra variables for profiles to 2 to store
817         ! rotation parameters
818         nveloextr = 2
819
820         jveloset = 0
821         
822         IF ( ln_velavcur ) jveloset = jveloset + 1
823         IF ( ln_velhrcur ) jveloset = jveloset + 1
824         IF ( ln_velavadcp ) jveloset = jveloset + 1
825         IF ( ln_velhradcp ) jveloset = jveloset + 1
826         IF (ln_velfb) jveloset = jveloset + jnumvelfb
827
828         nvelosets = jveloset
829         IF ( nvelosets > 0 ) THEN
830            ALLOCATE( velodata(nvelosets) )
831            ALLOCATE( veldatqc(nvelosets) )
832            ALLOCATE( ld_velav(nvelosets) )
833         ENDIF
834         
835         jveloset = 0
836         
837         ! Daily averaged data
838
839         IF ( ln_velavcur ) THEN
840           
841            jveloset = jveloset + 1
842           
843            ld_velav(jveloset) = .TRUE.
844           
845            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
846               &                  velavcurfiles(1:jnumvelavcur), &
847               &                  nvelovars, nveloextr, &
848               &                  nitend-nit000+2,              &
849               &                  dobsini, dobsend, ln_ignmis, &
850               &                  ld_velav(jveloset), &
851               &                  .FALSE. )
852           
853            DO jvar = 1, 2
854               CALL obs_prof_staend( velodata(jveloset), jvar )
855            END DO
856           
857            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
858               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
859           
860         ENDIF
861
862         ! High frequency data
863
864         IF ( ln_velhrcur ) THEN
865           
866            jveloset = jveloset + 1
867           
868            ld_velav(jveloset) = .FALSE.
869               
870            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
871               &                  velhrcurfiles(1:jnumvelhrcur), &
872               &                  nvelovars, nveloextr, &
873               &                  nitend-nit000+2,              &
874               &                  dobsini, dobsend, ln_ignmis, &
875               &                  ld_velav(jveloset), &
876               &                  .FALSE. )
877           
878            DO jvar = 1, 2
879               CALL obs_prof_staend( velodata(jveloset), jvar )
880            END DO
881           
882            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
883               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
884           
885         ENDIF
886
887         ! Daily averaged data
888
889         IF ( ln_velavadcp ) THEN
890           
891            jveloset = jveloset + 1
892           
893            ld_velav(jveloset) = .TRUE.
894           
895            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
896               &                  velavadcpfiles(1:jnumvelavadcp), &
897               &                  nvelovars, nveloextr, &
898               &                  nitend-nit000+2,              &
899               &                  dobsini, dobsend, ln_ignmis, &
900               &                  ld_velav(jveloset), &
901               &                  .FALSE. )
902           
903            DO jvar = 1, 2
904               CALL obs_prof_staend( velodata(jveloset), jvar )
905            END DO
906           
907            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
908               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
909           
910         ENDIF
911
912         ! High frequency data
913
914         IF ( ln_velhradcp ) THEN
915           
916            jveloset = jveloset + 1
917           
918            ld_velav(jveloset) = .FALSE.
919               
920            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
921               &                  velhradcpfiles(1:jnumvelhradcp), &
922               &                  nvelovars, nveloextr, &
923               &                  nitend-nit000+2,              &
924               &                  dobsini, dobsend, ln_ignmis, &
925               &                  ld_velav(jveloset), &
926               &                  .FALSE. )
927           
928            DO jvar = 1, 2
929               CALL obs_prof_staend( velodata(jveloset), jvar )
930            END DO
931           
932            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
933               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
934           
935         ENDIF
936
937         IF ( ln_velfb ) THEN
938
939            DO jset = 1, jnumvelfb
940           
941               jveloset = jveloset + 1
942
943               ld_velav(jveloset) = ln_velfb_av(jset)
944               
945               CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
946                  &                  velfbfiles(jset:jset), &
947                  &                  nvelovars, nveloextr, &
948                  &                  nitend-nit000+2,              &
949                  &                  dobsini, dobsend, ln_ignmis, &
950                  &                  ld_velav(jveloset), &
951                  &                  .FALSE. )
952               
953               DO jvar = 1, 2
954                  CALL obs_prof_staend( velodata(jveloset), jvar )
955               END DO
956               
957               CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
958                  &              ln_vel3d, ln_nea, ld_velav(jveloset) )
959
960
961            END DO
962           
963         ENDIF
964
965      ENDIF
966     
967   END SUBROUTINE dia_obs_init
968
969   SUBROUTINE dia_obs( kstp )
970      !!----------------------------------------------------------------------
971      !!                    ***  ROUTINE dia_obs  ***
972      !!         
973      !! ** Purpose : Call the observation operators on each time step
974      !!
975      !! ** Method  : Call the observation operators on each time step to
976      !!              compute the model equivalent of the following date:
977      !!               - T profiles
978      !!               - S profiles
979      !!               - Sea surface height (referenced to a mean)
980      !!               - Sea surface temperature
981      !!               - Sea surface salinity
982      !!               - Velocity component (U,V) profiles
983      !!
984      !! ** Action  :
985      !!
986      !! History :
987      !!        !  06-03  (K. Mogensen) Original code
988      !!        !  06-05  (K. Mogensen) Reformatted
989      !!        !  06-10  (A. Weaver) Cleaning
990      !!        !  07-03  (K. Mogensen) General handling of profiles
991      !!        !  07-04  (G. Smith) Generalized surface operators
992      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
993      !!----------------------------------------------------------------------
994      !! * Modules used
995      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
996         & rdt,           &                       
997         & gdept_1d,       &             
998         & tmask, umask, vmask                           
999      USE phycst, ONLY : &              ! Physical constants
1000         & rday                         
1001      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
1002         & tsn,  &             
1003         & un, vn,  &
1004         & sshn
1005#if defined  key_lim3
1006      USE ice, ONLY : &                     ! LIM Ice model variables
1007         & frld
1008#endif
1009#if defined key_lim2
1010      USE ice_2, ONLY : &                     ! LIM Ice model variables
1011         & frld
1012#endif
1013      IMPLICIT NONE
1014
1015      !! * Arguments
1016      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
1017      !! * Local declarations
1018      INTEGER :: idaystp                ! Number of timesteps per day
1019      INTEGER :: jprofset               ! Profile data set loop variable
1020      INTEGER :: jslaset                ! SLA data set loop variable
1021      INTEGER :: jsstset                ! SST data set loop variable
1022      INTEGER :: jseaiceset             ! sea ice data set loop variable
1023      INTEGER :: jveloset               ! velocity profile data loop variable
1024      INTEGER :: jvar                   ! Variable number   
1025#if ! defined key_lim2 && ! defined key_lim3
1026      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
1027#endif
1028      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1029 
1030#if ! defined key_lim2 && ! defined key_lim3
1031      CALL wrk_alloc(jpi,jpj,frld) 
1032#endif
1033
1034      IF(lwp) THEN
1035         WRITE(numout,*)
1036         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1037         WRITE(numout,*) '~~~~~~~'
1038      ENDIF
1039
1040      idaystp = NINT( rday / rdt )
1041
1042      !-----------------------------------------------------------------------
1043      ! No LIM => frld == 0.0_wp
1044      !-----------------------------------------------------------------------
1045#if ! defined key_lim2 && ! defined key_lim3
1046      frld(:,:) = 0.0_wp
1047#endif
1048      !-----------------------------------------------------------------------
1049      ! Depending on switches call various observation operators
1050      !-----------------------------------------------------------------------
1051
1052      !  - Temperature/salinity profiles
1053      IF ( ln_t3d .OR. ln_s3d ) THEN
1054         DO jprofset = 1, nprofsets
1055            IF ( ld_enact(jprofset) ) THEN
1056               CALL obs_pro_opt( prodatqc(jprofset),                     &
1057                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   &
1058                  &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &
1059                  &              gdept_1d, tmask, n1dint, n2dint,        &
1060                  &              kdailyavtypes = endailyavtypes )
1061            ELSE
1062               CALL obs_pro_opt( prodatqc(jprofset),                     &
1063                  &              kstp, jpi, jpj, jpk, nit000, idaystp,   &
1064                  &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   &
1065                  &              gdept_1d, tmask, n1dint, n2dint              )
1066            ENDIF
1067         END DO
1068      ENDIF
1069
1070      !  - Sea surface anomaly
1071      IF ( ln_sla ) THEN
1072         DO jslaset = 1, nslasets
1073            CALL obs_sla_opt( sladatqc(jslaset),            &
1074               &              kstp, jpi, jpj, nit000, sshn, &
1075               &              tmask(:,:,1), n2dint )
1076         END DO         
1077      ENDIF
1078
1079      !  - Sea surface temperature
1080      IF ( ln_sst ) THEN
1081         DO jsstset = 1, nsstsets
1082            CALL obs_sst_opt( sstdatqc(jsstset),                &
1083               &              kstp, jpi, jpj, nit000, idaystp,  &
1084               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  &
1085               &              n2dint, ld_sstnight(jsstset) )
1086         END DO
1087      ENDIF
1088
1089      !  - Sea surface salinity
1090      IF ( ln_sss ) THEN
1091         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1092      ENDIF
1093
1094#if defined key_lim2 || defined key_lim3
1095      IF ( ln_seaice ) THEN
1096         DO jseaiceset = 1, nseaicesets
1097            CALL obs_seaice_opt( seaicedatqc(jseaiceset),      &
1098               &              kstp, jpi, jpj, nit000, 1.-frld, &
1099               &              tmask(:,:,1), n2dint )
1100         END DO
1101      ENDIF     
1102#endif
1103
1104      !  - Velocity profiles
1105      IF ( ln_vel3d ) THEN
1106         DO jveloset = 1, nvelosets
1107           ! zonal component of velocity
1108           CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
1109              &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
1110                             n1dint, n2dint, ld_velav(jveloset) )
1111         END DO
1112      ENDIF
1113
1114#if ! defined key_lim2 && ! defined key_lim3
1115      CALL wrk_dealloc(jpi,jpj,frld) 
1116#endif
1117
1118   END SUBROUTINE dia_obs
1119 
1120   SUBROUTINE dia_obs_wri 
1121      !!----------------------------------------------------------------------
1122      !!                    ***  ROUTINE dia_obs_wri  ***
1123      !!         
1124      !! ** Purpose : Call observation diagnostic output routines
1125      !!
1126      !! ** Method  : Call observation diagnostic output routines
1127      !!
1128      !! ** Action  :
1129      !!
1130      !! History :
1131      !!        !  06-03  (K. Mogensen) Original code
1132      !!        !  06-05  (K. Mogensen) Reformatted
1133      !!        !  06-10  (A. Weaver) Cleaning
1134      !!        !  07-03  (K. Mogensen) General handling of profiles
1135      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
1136      !!----------------------------------------------------------------------
1137      IMPLICIT NONE
1138
1139      !! * Local declarations
1140
1141      INTEGER :: jprofset                 ! Profile data set loop variable
1142      INTEGER :: jveloset                 ! Velocity data set loop variable
1143      INTEGER :: jslaset                  ! SLA data set loop variable
1144      INTEGER :: jsstset                  ! SST data set loop variable
1145      INTEGER :: jseaiceset               ! Sea Ice data set loop variable
1146      INTEGER :: jset
1147      INTEGER :: jfbini
1148      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1149      CHARACTER(LEN=10) :: cdtmp
1150      !-----------------------------------------------------------------------
1151      ! Depending on switches call various observation output routines
1152      !-----------------------------------------------------------------------
1153
1154      !  - Temperature/salinity profiles
1155
1156      IF( ln_t3d .OR. ln_s3d ) THEN
1157
1158         ! Copy data from prodatqc to profdata structures
1159         DO jprofset = 1, nprofsets
1160
1161            CALL obs_prof_decompress( prodatqc(jprofset), &
1162                 &                    profdata(jprofset), .TRUE., numout )
1163
1164         END DO
1165
1166         ! Write the profiles.
1167
1168         jprofset = 0
1169
1170         ! ENACT insitu data
1171
1172         IF ( ln_ena ) THEN
1173           
1174            jprofset = jprofset + 1
1175
1176            CALL obs_wri_p3d( 'enact', profdata(jprofset) )
1177
1178         ENDIF
1179
1180         ! Coriolis insitu data
1181
1182         IF ( ln_cor ) THEN
1183           
1184            jprofset = jprofset + 1
1185
1186            CALL obs_wri_p3d( 'corio', profdata(jprofset) )
1187           
1188         ENDIF
1189         
1190         ! Feedback insitu data
1191
1192         IF ( ln_profb ) THEN
1193
1194            jfbini = jprofset + 1
1195
1196            DO jprofset = jfbini, nprofsets
1197               
1198               jset = jprofset - jfbini + 1
1199               WRITE(cdtmp,'(A,I2.2)')'profb_',jset
1200               CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
1201
1202            END DO
1203
1204         ENDIF
1205
1206      ENDIF
1207
1208      !  - Sea surface anomaly
1209      IF ( ln_sla ) THEN
1210
1211         ! Copy data from sladatqc to sladata structures
1212         DO jslaset = 1, nslasets
1213
1214              CALL obs_surf_decompress( sladatqc(jslaset), &
1215                 &                    sladata(jslaset), .TRUE., numout )
1216
1217         END DO
1218
1219         jslaset = 0 
1220
1221         ! Write the AVISO SLA data
1222
1223         IF ( ln_sladt ) THEN
1224           
1225            jslaset = 1
1226            CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
1227            jslaset = 2
1228            CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
1229
1230         ENDIF
1231
1232         IF ( ln_slafb ) THEN
1233           
1234            jfbini = jslaset + 1
1235
1236            DO jslaset = jfbini, nslasets
1237               
1238               jset = jslaset - jfbini + 1
1239               WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
1240               CALL obs_wri_sla( cdtmp, sladata(jslaset) )
1241
1242            END DO
1243
1244         ENDIF
1245
1246      ENDIF
1247
1248      !  - Sea surface temperature
1249      IF ( ln_sst ) THEN
1250
1251         ! Copy data from sstdatqc to sstdata structures
1252         DO jsstset = 1, nsstsets
1253     
1254              CALL obs_surf_decompress( sstdatqc(jsstset), &
1255                 &                    sstdata(jsstset), .TRUE., numout )
1256
1257         END DO
1258
1259         jsstset = 0 
1260
1261         ! Write the AVISO SST data
1262
1263         IF ( ln_reysst ) THEN
1264           
1265            jsstset = jsstset + 1
1266            CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
1267
1268         ENDIF
1269
1270         IF ( ln_ghrsst ) THEN
1271           
1272            jsstset = jsstset + 1
1273            CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
1274
1275         ENDIF
1276
1277         IF ( ln_sstfb ) THEN
1278           
1279            jfbini = jsstset + 1
1280
1281            DO jsstset = jfbini, nsstsets
1282               
1283               jset = jsstset - jfbini + 1
1284               WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
1285               CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
1286
1287            END DO
1288
1289         ENDIF
1290
1291      ENDIF
1292
1293      !  - Sea surface salinity
1294      IF ( ln_sss ) THEN
1295         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1296      ENDIF
1297
1298      !  - Sea Ice Concentration
1299      IF ( ln_seaice ) THEN
1300
1301         ! Copy data from seaicedatqc to seaicedata structures
1302         DO jseaiceset = 1, nseaicesets
1303
1304              CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
1305                 &                    seaicedata(jseaiceset), .TRUE., numout )
1306
1307         END DO
1308
1309         ! Write the Sea Ice data
1310         DO jseaiceset = 1, nseaicesets
1311     
1312            WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
1313            CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
1314
1315         END DO
1316
1317      ENDIF
1318     
1319      ! Velocity data
1320      IF( ln_vel3d ) THEN
1321
1322         ! Copy data from veldatqc to velodata structures
1323         DO jveloset = 1, nvelosets
1324
1325            CALL obs_prof_decompress( veldatqc(jveloset), &
1326                 &                    velodata(jveloset), .TRUE., numout )
1327
1328         END DO
1329
1330         ! Write the profiles.
1331
1332         jveloset = 0
1333
1334         ! Daily averaged data
1335
1336         IF ( ln_velavcur ) THEN
1337           
1338            jveloset = jveloset + 1
1339
1340            CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
1341
1342         ENDIF
1343
1344         ! High frequency data
1345
1346         IF ( ln_velhrcur ) THEN
1347           
1348            jveloset = jveloset + 1
1349
1350            CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
1351
1352         ENDIF
1353
1354         ! Daily averaged data
1355
1356         IF ( ln_velavadcp ) THEN
1357           
1358            jveloset = jveloset + 1
1359
1360            CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
1361
1362         ENDIF
1363
1364         ! High frequency data
1365
1366         IF ( ln_velhradcp ) THEN
1367           
1368            jveloset = jveloset + 1
1369           
1370            CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
1371               
1372         ENDIF
1373
1374         ! Feedback velocity data
1375
1376         IF ( ln_velfb ) THEN
1377
1378            jfbini = jveloset + 1
1379
1380            DO jveloset = jfbini, nvelosets
1381               
1382               jset = jveloset - jfbini + 1
1383               WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
1384               CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
1385
1386            END DO
1387
1388         ENDIF
1389         
1390      ENDIF
1391
1392   END SUBROUTINE dia_obs_wri
1393
1394   SUBROUTINE dia_obs_dealloc
1395      IMPLICIT NONE
1396      !!----------------------------------------------------------------------
1397      !!                    *** ROUTINE dia_obs_dealloc ***
1398      !!
1399      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1400      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1401      !!
1402      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1403      !!
1404      !!  ** Action :
1405      !!
1406      !!----------------------------------------------------------------------
1407      !! obs_grid deallocation
1408      CALL obs_grid_deallocate
1409
1410      !! diaobs deallocation
1411      IF ( nprofsets > 0 ) THEN
1412          DEALLOCATE(ld_enact, &
1413                  &  profdata, &
1414                  &  prodatqc)
1415      END IF
1416      IF ( ln_sla ) THEN
1417          DEALLOCATE(sladata, &
1418                  &  sladatqc)
1419      END IF
1420      IF ( ln_seaice ) THEN
1421          DEALLOCATE(sladata, &
1422                  &  sladatqc)
1423      END IF
1424      IF ( ln_sst ) THEN
1425          DEALLOCATE(sstdata, &
1426                  &  sstdatqc)
1427      END IF
1428      IF ( ln_vel3d ) THEN
1429          DEALLOCATE(ld_velav, &
1430                  &  velodata, &
1431                  &  veldatqc)
1432      END IF
1433   END SUBROUTINE dia_obs_dealloc
1434
1435   SUBROUTINE ini_date( ddobsini )
1436      !!----------------------------------------------------------------------
1437      !!                    ***  ROUTINE ini_date  ***
1438      !!         
1439      !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
1440      !!
1441      !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format
1442      !!
1443      !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format
1444      !!
1445      !! History :
1446      !!        !  06-03  (K. Mogensen)  Original code
1447      !!        !  06-05  (K. Mogensen)  Reformatted
1448      !!        !  06-10  (A. Weaver) Cleaning
1449      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1450      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1451      !!----------------------------------------------------------------------
1452      USE phycst, ONLY : &            ! Physical constants
1453         & rday
1454!      USE daymod, ONLY : &            ! Time variables
1455!         & nmonth_len           
1456      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1457         & rdt
1458
1459      IMPLICIT NONE
1460
1461      !! * Arguments
1462      REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS
1463
1464      !! * Local declarations
1465      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1466      INTEGER :: imon
1467      INTEGER :: iday
1468      INTEGER :: ihou
1469      INTEGER :: imin
1470      INTEGER :: imday         ! Number of days in month.
1471      REAL(KIND=wp) :: zdayfrc ! Fraction of day
1472
1473      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1474
1475      !!----------------------------------------------------------------------
1476      !! Initial date initialization (year, month, day, hour, minute)
1477      !! (This assumes that the initial date is for 00z))
1478      !!----------------------------------------------------------------------
1479      iyea =   ndate0 / 10000
1480      imon = ( ndate0 - iyea * 10000 ) / 100
1481      iday =   ndate0 - iyea * 10000 - imon * 100
1482      ihou = 0
1483      imin = 0
1484
1485      !!----------------------------------------------------------------------
1486      !! Compute number of days + number of hours + min since initial time
1487      !!----------------------------------------------------------------------
1488      iday = iday + ( nit000 -1 ) * rdt / rday
1489      zdayfrc = ( nit000 -1 ) * rdt / rday
1490      zdayfrc = zdayfrc - aint(zdayfrc)
1491      ihou = int( zdayfrc * 24 )
1492      imin = int( (zdayfrc * 24 - ihou) * 60 )
1493
1494      !!-----------------------------------------------------------------------
1495      !! Convert number of days (iday) into a real date
1496      !!----------------------------------------------------------------------
1497
1498      CALL calc_month_len( iyea, imonth_len )
1499     
1500      DO WHILE ( iday > imonth_len(imon) )
1501         iday = iday - imonth_len(imon)
1502         imon = imon + 1 
1503         IF ( imon > 12 ) THEN
1504            imon = 1
1505            iyea = iyea + 1
1506            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1507         ENDIF
1508      END DO
1509
1510      !!----------------------------------------------------------------------
1511      !! Convert it into YYYYMMDD.HHMMSS format.
1512      !!----------------------------------------------------------------------
1513      ddobsini = iyea * 10000_dp + imon * 100_dp + &
1514         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
1515
1516
1517   END SUBROUTINE ini_date
1518
1519   SUBROUTINE fin_date( ddobsfin )
1520      !!----------------------------------------------------------------------
1521      !!                    ***  ROUTINE fin_date  ***
1522      !!         
1523      !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
1524      !!
1525      !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format
1526      !!
1527      !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format
1528      !!
1529      !! History :
1530      !!        !  06-03  (K. Mogensen)  Original code
1531      !!        !  06-05  (K. Mogensen)  Reformatted
1532      !!        !  06-10  (A. Weaver) Cleaning
1533      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1534      !!----------------------------------------------------------------------
1535      USE phycst, ONLY : &            ! Physical constants
1536         & rday
1537!      USE daymod, ONLY : &            ! Time variables
1538!         & nmonth_len               
1539      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1540         & rdt
1541
1542      IMPLICIT NONE
1543
1544      !! * Arguments
1545      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
1546
1547      !! * Local declarations
1548      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1549      INTEGER :: imon
1550      INTEGER :: iday
1551      INTEGER :: ihou
1552      INTEGER :: imin
1553      INTEGER :: imday         ! Number of days in month.
1554      REAL(KIND=wp) :: zdayfrc       ! Fraction of day
1555         
1556      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1557           
1558      !-----------------------------------------------------------------------
1559      ! Initial date initialization (year, month, day, hour, minute)
1560      ! (This assumes that the initial date is for 00z)
1561      !-----------------------------------------------------------------------
1562      iyea =   ndate0 / 10000
1563      imon = ( ndate0 - iyea * 10000 ) / 100
1564      iday =   ndate0 - iyea * 10000 - imon * 100
1565      ihou = 0
1566      imin = 0
1567     
1568      !-----------------------------------------------------------------------
1569      ! Compute number of days + number of hours + min since initial time
1570      !-----------------------------------------------------------------------
1571      iday    = iday +  nitend  * rdt / rday
1572      zdayfrc =  nitend  * rdt / rday
1573      zdayfrc = zdayfrc - AINT( zdayfrc )
1574      ihou    = INT( zdayfrc * 24 )
1575      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
1576
1577      !-----------------------------------------------------------------------
1578      ! Convert number of days (iday) into a real date
1579      !----------------------------------------------------------------------
1580
1581      CALL calc_month_len( iyea, imonth_len )
1582     
1583      DO WHILE ( iday > imonth_len(imon) )
1584         iday = iday - imonth_len(imon)
1585         imon = imon + 1 
1586         IF ( imon > 12 ) THEN
1587            imon = 1
1588            iyea = iyea + 1
1589            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1590         ENDIF
1591      END DO
1592
1593      !-----------------------------------------------------------------------
1594      ! Convert it into YYYYMMDD.HHMMSS format
1595      !-----------------------------------------------------------------------
1596      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
1597         &     + ihou * 0.01_dp  + imin * 0.0001_dp
1598
1599    END SUBROUTINE fin_date
1600   
1601END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.