New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diaobs.F90 in branches/2014/dev_r4650_UKMO14.11_SETTE_OBSASM/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2014/dev_r4650_UKMO14.11_SETTE_OBSASM/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 4751

Last change on this file since 4751 was 4751, checked in by djlea, 10 years ago

Changes to include an OBS test in SETTE. At the moment this uses an example profile observation.

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