source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6016

Last change on this file since 6016 was 6016, checked in by kingr, 5 years ago

#1642 bug fixes for general vertical coord obsoper

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