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/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_surf_bgc/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6855

Last change on this file since 6855 was 6855, checked in by dford, 8 years ago

Initial implementation of observation operator for SPM.

File size: 74.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_const, ONLY: obfillflt ! Fill value
24   USE obs_fbm, ONLY: ln_cl4    ! Class 4 diagnostic switch
25   USE obs_read_prof            ! Reading and allocation of observations (Coriolis)
26   USE obs_read_sla             ! Reading and allocation of SLA observations 
27   USE obs_read_sst             ! Reading and allocation of SST observations 
28   USE obs_sstbias              ! Bias correction routine for SST
29   USE obs_readmdt              ! Reading and allocation of MDT for SLA.
30   USE obs_read_seaice          ! Reading and allocation of Sea Ice observations 
31   USE obs_read_vel             ! Reading and allocation of velocity component observations
32   USE obs_read_logchl          ! Reading and allocation of logchl observations
33   USE obs_read_spm             ! Reading and allocation of spm observations
34   USE obs_prep                 ! Preparation of obs. (grid search etc).
35   USE obs_oper                 ! Observation operators
36   USE obs_write                ! Writing of observation related diagnostics
37   USE obs_grid                 ! Grid searching
38   USE obs_read_altbias         ! Bias treatment for altimeter
39   USE obs_profiles_def         ! Profile data definitions
40   USE obs_profiles             ! Profile data storage
41   USE obs_surf_def             ! Surface data definitions
42   USE obs_sla                  ! SLA data storage
43   USE obs_sst                  ! SST data storage
44   USE obs_seaice               ! Sea Ice data storage
45   USE obs_logchl               ! logchl data storage
46   USE obs_spm                  ! spm data storage
47   USE obs_types                ! Definitions for observation types
48   USE mpp_map                  ! MPP mapping
49   USE lib_mpp                  ! For ctl_warn/stop
50
51   IMPLICIT NONE
52
53   !! * Routine accessibility
54   PRIVATE
55   PUBLIC dia_obs_init, &  ! Initialize and read observations
56      &   dia_obs,      &  ! Compute model equivalent to observations
57      &   dia_obs_wri,  &  ! Write model equivalent to observations
58      &   dia_obs_dealloc  ! Deallocate dia_obs data
59
60   !! * Shared Module variables
61   LOGICAL, PUBLIC, PARAMETER :: &
62#if defined key_diaobs
63      & lk_diaobs = .TRUE.   !: Logical switch for observation diangostics
64#else
65      & lk_diaobs = .FALSE.  !: Logical switch for observation diangostics
66#endif
67
68   !! * Module variables
69   LOGICAL, PUBLIC :: ln_t3d         !: Logical switch for temperature profiles
70   LOGICAL, PUBLIC :: ln_s3d         !: Logical switch for salinity profiles
71   LOGICAL, PUBLIC :: ln_ena         !: Logical switch for the ENACT data set
72   LOGICAL, PUBLIC :: ln_cor         !: Logical switch for the Coriolis data set
73   LOGICAL, PUBLIC :: ln_profb       !: Logical switch for profile feedback datafiles
74   LOGICAL, PUBLIC :: ln_sla         !: Logical switch for sea level anomalies
75   LOGICAL, PUBLIC :: ln_sladt       !: Logical switch for SLA from AVISO files
76   LOGICAL, PUBLIC :: ln_slafb       !: Logical switch for SLA from feedback files
77   LOGICAL, PUBLIC :: ln_sst         !: Logical switch for sea surface temperature
78   LOGICAL, PUBLIC :: ln_reysst      !: Logical switch for Reynolds sea surface temperature
79   LOGICAL, PUBLIC :: ln_ghrsst      !: Logical switch for GHRSST data
80   LOGICAL, PUBLIC :: ln_sstfb       !: Logical switch for SST from feedback files
81   LOGICAL, PUBLIC :: ln_seaice      !: Logical switch for sea ice concentration
82   LOGICAL, PUBLIC :: ln_vel3d       !: Logical switch for velocity component (u,v) observations
83   LOGICAL, PUBLIC :: ln_velavcur    !: Logical switch for raw daily averaged netCDF current meter vel. data
84   LOGICAL, PUBLIC :: ln_velhrcur    !: Logical switch for raw high freq netCDF current meter vel. data
85   LOGICAL, PUBLIC :: ln_velavadcp   !: Logical switch for raw daily averaged netCDF ADCP vel. data
86   LOGICAL, PUBLIC :: ln_velhradcp   !: Logical switch for raw high freq netCDF ADCP vel. data
87   LOGICAL, PUBLIC :: ln_velfb       !: Logical switch for velocities from feedback files
88   LOGICAL, PUBLIC :: ln_logchl      !: Logical switch for log10(chlorophyll)
89   LOGICAL, PUBLIC :: ln_logchlfb    !: Logical switch for logchl from feedback files
90   LOGICAL, PUBLIC :: ln_spm         !: Logical switch for spm
91   LOGICAL, PUBLIC :: ln_spmfb       !: Logical switch for spm from feedback files
92   LOGICAL, PUBLIC :: ln_ssh         !: Logical switch for sea surface height
93   LOGICAL, PUBLIC :: ln_sss         !: Logical switch for sea surface salinity
94   LOGICAL, PUBLIC :: ln_sstnight    !: Logical switch for night mean SST observations
95   LOGICAL, PUBLIC :: ln_nea         !: Remove observations near land
96   LOGICAL, PUBLIC :: ln_altbias     !: Logical switch for altimeter bias 
97   LOGICAL, PUBLIC :: ln_ignmis      !: Logical switch for ignoring missing files
98   LOGICAL, PUBLIC :: ln_s_at_t      !: Logical switch to compute model S at T observations
99   LOGICAL, PUBLIC :: ln_sstbias     !: Logical switch for bias corection of SST
100
101   REAL(KIND=dp), PUBLIC :: dobsini   !: Observation window start date YYYYMMDD.HHMMSS
102   REAL(KIND=dp), PUBLIC :: dobsend   !: Observation window end date YYYYMMDD.HHMMSS
103 
104   INTEGER, PUBLIC :: n1dint       !: Vertical interpolation method
105   INTEGER, PUBLIC :: n2dint       !: Horizontal interpolation method
106
107   INTEGER, DIMENSION(imaxavtypes) :: &
108      & endailyavtypes !: ENACT data types which are daily average
109
110   INTEGER, PARAMETER :: MaxNumFiles = 1000
111   LOGICAL, DIMENSION(MaxNumFiles) :: &
112      & ln_profb_ena, & !: Is the feedback files from ENACT data ?
113   !                    !: If so use endailyavtypes
114      & ln_profb_enatim !: Change tim for 820 enact data set.
115   
116   INTEGER, DIMENSION(MaxNumFiles), PUBLIC :: sstbias_type !SST bias type
117
118   LOGICAL, DIMENSION(MaxNumFiles) :: &
119      & ln_velfb_av   !: Is the velocity feedback files daily average?
120   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
121      & ld_enact     !: Profile data is ENACT so use endailyavtypes
122   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
123      & ld_velav     !: Velocity data is daily averaged
124   LOGICAL, DIMENSION(:), ALLOCATABLE :: &
125      & ld_sstnight  !: SST observation corresponds to night mean
126
127   !!----------------------------------------------------------------------
128   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
129   !! $Id$
130   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
131   !!----------------------------------------------------------------------
132
133   !! * Substitutions
134#  include "domzgr_substitute.h90"
135CONTAINS
136
137   SUBROUTINE dia_obs_init
138      !!----------------------------------------------------------------------
139      !!                    ***  ROUTINE dia_obs_init  ***
140      !!         
141      !! ** Purpose : Initialize and read observations
142      !!
143      !! ** Method  : Read the namelist and call reading routines
144      !!
145      !! ** Action  : Read the namelist and call reading routines
146      !!
147      !! History :
148      !!        !  06-03  (K. Mogensen) Original code
149      !!        !  06-05  (A. Weaver) Reformatted
150      !!        !  06-10  (A. Weaver) Cleaning and add controls
151      !!        !  07-03  (K. Mogensen) General handling of profiles
152      !!        !  14-08  (J.While) Incorporated SST bias correction
153      !!----------------------------------------------------------------------
154
155      IMPLICIT NONE
156
157      !! * Local declarations
158      CHARACTER(len=128) :: enactfiles(MaxNumFiles)
159      CHARACTER(len=128) :: coriofiles(MaxNumFiles)
160      CHARACTER(len=128) :: profbfiles(MaxNumFiles)
161      CHARACTER(len=128) :: sstfiles(MaxNumFiles)     
162      CHARACTER(len=128) :: sstfbfiles(MaxNumFiles)
163      CHARACTER(len=128) :: sstbias_files(MaxNumFiles) 
164      CHARACTER(len=128) :: slafilesact(MaxNumFiles)     
165      CHARACTER(len=128) :: slafilespas(MaxNumFiles)     
166      CHARACTER(len=128) :: slafbfiles(MaxNumFiles)
167      CHARACTER(len=128) :: seaicefiles(MaxNumFiles)           
168      CHARACTER(len=128) :: velcurfiles(MaxNumFiles) 
169      CHARACTER(len=128) :: veladcpfiles(MaxNumFiles)   
170      CHARACTER(len=128) :: velavcurfiles(MaxNumFiles)
171      CHARACTER(len=128) :: velhrcurfiles(MaxNumFiles)
172      CHARACTER(len=128) :: velavadcpfiles(MaxNumFiles)
173      CHARACTER(len=128) :: velhradcpfiles(MaxNumFiles)
174      CHARACTER(len=128) :: velfbfiles(MaxNumFiles)
175      CHARACTER(len=128) :: logchlfiles(MaxNumFiles)
176      CHARACTER(len=128) :: logchlfbfiles(MaxNumFiles)
177      CHARACTER(len=128) :: spmfiles(MaxNumFiles)
178      CHARACTER(len=128) :: spmfbfiles(MaxNumFiles)
179      CHARACTER(LEN=128) :: reysstname
180      CHARACTER(LEN=12)  :: reysstfmt
181      CHARACTER(LEN=128) :: bias_file
182      CHARACTER(LEN=20)  :: datestr=" ", timestr=" "
183      NAMELIST/namobs/ln_ena, ln_cor, ln_profb, ln_t3d, ln_s3d,       &
184         &            ln_sla, ln_sladt, ln_slafb,                     &
185         &            ln_ssh, ln_sst, ln_sstfb, ln_sss, ln_nea,       &
186         &            enactfiles, coriofiles, profbfiles,             &
187         &            slafilesact, slafilespas, slafbfiles,           &
188         &            sstfiles, sstfbfiles,                           &
189         &            ln_seaice, seaicefiles,                         &
190         &            dobsini, dobsend, n1dint, n2dint,               &
191         &            nmsshc, mdtcorr, mdtcutoff,                     &
192         &            ln_reysst, ln_ghrsst, reysstname, reysstfmt,    &
193         &            ln_sstnight,                                    &
194         &            ln_grid_search_lookup,                          &
195         &            grid_search_file, grid_search_res,              &
196         &            ln_grid_global, bias_file, ln_altbias,          &
197         &            endailyavtypes, ln_s_at_t, ln_profb_ena,        &
198         &            ln_vel3d, ln_velavcur, velavcurfiles,           &
199         &            ln_velhrcur, velhrcurfiles,                     &
200         &            ln_velavadcp, velavadcpfiles,                   &
201         &            ln_velhradcp, velhradcpfiles,                   &
202         &            ln_velfb, velfbfiles, ln_velfb_av,              &
203         &            ln_logchl, ln_logchlfb,                         &
204         &            logchlfiles, logchlfbfiles,                     &
205         &            ln_spm, ln_spmfb,                               &
206         &            spmfiles, spmfbfiles,                           &
207         &            ln_profb_enatim, ln_ignmis, ln_cl4,             &
208         &            ln_sstbias, sstbias_files
209
210      INTEGER :: jprofset
211      INTEGER :: jveloset
212      INTEGER :: jvar
213      INTEGER :: jnumenact
214      INTEGER :: jnumcorio
215      INTEGER :: jnumprofb
216      INTEGER :: jnumslaact
217      INTEGER :: jnumslapas
218      INTEGER :: jnumslafb
219      INTEGER :: jnumsst
220      INTEGER :: jnumsstfb
221      INTEGER :: jnumsstbias
222      INTEGER :: jnumseaice
223      INTEGER :: jnumvelavcur
224      INTEGER :: jnumvelhrcur 
225      INTEGER :: jnumvelavadcp
226      INTEGER :: jnumvelhradcp   
227      INTEGER :: jnumvelfb
228      INTEGER :: jnumlogchl
229      INTEGER :: jnumlogchlfb
230      INTEGER :: jnumspm
231      INTEGER :: jnumspmfb
232      INTEGER :: ji
233      INTEGER :: jset
234      INTEGER :: ios                 ! Local integer output status for namelist read
235      LOGICAL :: lmask(MaxNumFiles), ll_u3d, ll_v3d
236
237      !-----------------------------------------------------------------------
238      ! Read namelist parameters
239      !-----------------------------------------------------------------------
240
241      ln_logchl   = .FALSE.
242      ln_logchlfb = .FALSE.
243      ln_spm      = .FALSE.
244      ln_spmfb    = .FALSE.
245     
246      !Initalise all values in namelist arrays
247      enactfiles(:) = ''
248      coriofiles(:) = ''
249      profbfiles(:) = ''
250      slafilesact(:) = ''
251      slafilespas(:) = ''
252      slafbfiles(:) = ''
253      sstfiles(:)   = ''
254      sstfbfiles(:) = ''
255      seaicefiles(:) = ''
256      velcurfiles(:) = ''
257      veladcpfiles(:) = ''
258      velavcurfiles(:) = ''
259      velhrcurfiles(:) = ''
260      velavadcpfiles(:) = ''
261      velhradcpfiles(:) = ''
262      velfbfiles(:) = ''
263      velcurfiles(:) = ''
264      veladcpfiles(:) = ''
265      logchlfiles(:) = ''
266      logchlfbfiles(:) = ''
267      spmfiles(:) = ''
268      spmfbfiles(:) = ''
269      sstbias_files(:) = ''
270      endailyavtypes(:) = -1
271      endailyavtypes(1) = 820
272      ln_profb_ena(:) = .FALSE.
273      ln_profb_enatim(:) = .TRUE.
274      ln_velfb_av(:) = .FALSE.
275      ln_ignmis = .FALSE.
276
277      ! Read Namelist namobs : control observation diagnostics
278      REWIND( numnam_ref )              ! Namelist namobs in reference namelist : Diagnostic: control observation
279      READ  ( numnam_ref, namobs, IOSTAT = ios, ERR = 901)
280901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in reference namelist', lwp )
281
282      CALL ini_date( dobsini )
283      CALL fin_date( dobsend )
284 
285      REWIND( numnam_cfg )              ! Namelist namobs in configuration namelist : Diagnostic: control observation
286      READ  ( numnam_cfg, namobs, IOSTAT = ios, ERR = 902 )
287902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namobs in configuration namelist', lwp )
288      IF(lwm) WRITE ( numond, namobs )
289
290      ! Count number of files for each type
291      IF (ln_ena) THEN
292         lmask(:) = .FALSE.
293         WHERE (enactfiles(:) /= '') lmask(:) = .TRUE.
294         jnumenact = COUNT(lmask)
295      ENDIF
296      IF (ln_cor) THEN
297         lmask(:) = .FALSE.
298         WHERE (coriofiles(:) /= '') lmask(:) = .TRUE.
299         jnumcorio = COUNT(lmask)
300      ENDIF
301      IF (ln_profb) THEN
302         lmask(:) = .FALSE.
303         WHERE (profbfiles(:) /= '') lmask(:) = .TRUE.
304         jnumprofb = COUNT(lmask)
305      ENDIF
306      IF (ln_sladt) THEN
307         lmask(:) = .FALSE.
308         WHERE (slafilesact(:) /= '') lmask(:) = .TRUE.
309         jnumslaact = COUNT(lmask)
310         lmask(:) = .FALSE.
311         WHERE (slafilespas(:) /= '') lmask(:) = .TRUE.
312         jnumslapas = COUNT(lmask)
313      ENDIF
314      IF (ln_slafb) THEN
315         lmask(:) = .FALSE.
316         WHERE (slafbfiles(:) /= '') lmask(:) = .TRUE.
317         jnumslafb = COUNT(lmask)
318         lmask(:) = .FALSE.
319      ENDIF
320      IF (ln_ghrsst) THEN
321         lmask(:) = .FALSE.
322         WHERE (sstfiles(:) /= '') lmask(:) = .TRUE.
323         jnumsst = COUNT(lmask)
324      ENDIF     
325      IF (ln_sstfb) THEN
326         lmask(:) = .FALSE.
327         WHERE (sstfbfiles(:) /= '') lmask(:) = .TRUE.
328         jnumsstfb = COUNT(lmask)
329         lmask(:) = .FALSE.
330      ENDIF
331      IF (ln_sstbias) THEN 
332         lmask(:) = .FALSE. 
333         WHERE (sstbias_files(:) /= '') lmask(:) = .TRUE. 
334         jnumsstbias = COUNT(lmask) 
335         lmask(:) = .FALSE. 
336      ENDIF       
337      IF (ln_seaice) THEN
338         lmask(:) = .FALSE.
339         WHERE (seaicefiles(:) /= '') lmask(:) = .TRUE.
340         jnumseaice = COUNT(lmask)
341      ENDIF
342      IF (ln_velavcur) THEN
343         lmask(:) = .FALSE.
344         WHERE (velavcurfiles(:) /= '') lmask(:) = .TRUE.
345         jnumvelavcur = COUNT(lmask)
346      ENDIF
347      IF (ln_velhrcur) THEN
348         lmask(:) = .FALSE.
349         WHERE (velhrcurfiles(:) /= '') lmask(:) = .TRUE.
350         jnumvelhrcur = COUNT(lmask)
351      ENDIF
352      IF (ln_velavadcp) THEN
353         lmask(:) = .FALSE.
354         WHERE (velavadcpfiles(:) /= '') lmask(:) = .TRUE.
355         jnumvelavadcp = COUNT(lmask)
356      ENDIF
357      IF (ln_velhradcp) THEN
358         lmask(:) = .FALSE.
359         WHERE (velhradcpfiles(:) /= '') lmask(:) = .TRUE.
360         jnumvelhradcp = COUNT(lmask)
361      ENDIF
362      IF (ln_velfb) THEN
363         lmask(:) = .FALSE.
364         WHERE (velfbfiles(:) /= '') lmask(:) = .TRUE.
365         jnumvelfb = COUNT(lmask)
366         lmask(:) = .FALSE.
367      ENDIF
368      IF (ln_logchl) THEN
369         lmask(:) = .FALSE.
370         WHERE (logchlfiles(:) /= '') lmask(:) = .TRUE.
371         jnumlogchl = COUNT(lmask)
372      ENDIF
373      IF (ln_logchlfb) THEN
374         lmask(:) = .FALSE.
375         WHERE (logchlfbfiles(:) /= '') lmask(:) = .TRUE.
376         jnumlogchlfb = COUNT(lmask)
377      ENDIF
378      IF (ln_spm) THEN
379         lmask(:) = .FALSE.
380         WHERE (spmfiles(:) /= '') lmask(:) = .TRUE.
381         jnumspm = COUNT(lmask)
382      ENDIF
383      IF (ln_spmfb) THEN
384         lmask(:) = .FALSE.
385         WHERE (spmfbfiles(:) /= '') lmask(:) = .TRUE.
386         jnumspmfb = COUNT(lmask)
387      ENDIF
388     
389      ! Control print
390      IF(lwp) THEN
391         WRITE(numout,*)
392         WRITE(numout,*) 'dia_obs_init : Observation diagnostic initialization'
393         WRITE(numout,*) '~~~~~~~~~~~~'
394         WRITE(numout,*) '          Namelist namobs : set observation diagnostic parameters' 
395         WRITE(numout,*) '             Logical switch for T profile observations          ln_t3d = ', ln_t3d
396         WRITE(numout,*) '             Logical switch for S profile observations          ln_s3d = ', ln_s3d
397         WRITE(numout,*) '             Logical switch for ENACT insitu data set           ln_ena = ', ln_ena
398         WRITE(numout,*) '             Logical switch for Coriolis insitu data set        ln_cor = ', ln_cor
399         WRITE(numout,*) '             Logical switch for feedback insitu data set      ln_profb = ', ln_profb
400         WRITE(numout,*) '             Logical switch for SLA observations                ln_sla = ', ln_sla
401         WRITE(numout,*) '             Logical switch for AVISO SLA data                ln_sladt = ', ln_sladt
402         WRITE(numout,*) '             Logical switch for feedback SLA data             ln_slafb = ', ln_slafb
403         WRITE(numout,*) '             Logical switch for SSH observations                ln_ssh = ', ln_ssh
404         WRITE(numout,*) '             Logical switch for SST observations                ln_sst = ', ln_sst
405         WRITE(numout,*) '             Logical switch for Reynolds observations        ln_reysst = ', ln_reysst   
406         WRITE(numout,*) '             Logical switch for GHRSST observations          ln_ghrsst = ', ln_ghrsst
407         WRITE(numout,*) '             Logical switch for feedback SST data             ln_sstfb = ', ln_sstfb
408         WRITE(numout,*) '             Logical switch for SST bias correction         ln_sstbias = ', ln_sstbias
409         WRITE(numout,*) '             Logical switch for night-time SST obs         ln_sstnight = ', ln_sstnight
410         WRITE(numout,*) '             Logical switch for SSS observations                ln_sss = ', ln_sss
411         WRITE(numout,*) '             Logical switch for Sea Ice observations         ln_seaice = ', ln_seaice
412         WRITE(numout,*) '             Logical switch for velocity observations         ln_vel3d = ', ln_vel3d
413         WRITE(numout,*) '             Logical switch for velocity daily av. cur.    ln_velavcur = ', ln_velavcur
414         WRITE(numout,*) '             Logical switch for velocity high freq. cur.   ln_velhrcur = ', ln_velhrcur
415         WRITE(numout,*) '             Logical switch for velocity daily av. ADCP   ln_velavadcp = ', ln_velavadcp
416         WRITE(numout,*) '             Logical switch for velocity high freq. ADCP  ln_velhradcp = ', ln_velhradcp
417         WRITE(numout,*) '             Logical switch for feedback velocity data        ln_velfb = ', ln_velfb
418         WRITE(numout,*) '             Logical switch for logchl observations          ln_logchl = ', ln_logchl
419         WRITE(numout,*) '             Logical switch for feedback logchl data       ln_logchlfb = ', ln_logchlfb
420         WRITE(numout,*) '             Logical switch for spm observations                ln_spm = ', ln_spm
421         WRITE(numout,*) '             Logical switch for feedback spm data             ln_spmfb = ', ln_spmfb
422         WRITE(numout,*) '             Global distribtion of observations         ln_grid_global = ',ln_grid_global
423         WRITE(numout,*) &
424   '             Logical switch for obs grid search w/lookup table  ln_grid_search_lookup = ',ln_grid_search_lookup
425         IF (ln_grid_search_lookup) &
426            WRITE(numout,*) '             Grid search lookup file header       grid_search_file = ', grid_search_file
427         IF (ln_ena) THEN
428            DO ji = 1, jnumenact
429               WRITE(numout,'(1X,2A)') '             ENACT input observation file name          enactfiles = ', &
430                  TRIM(enactfiles(ji))
431            END DO
432         ENDIF
433         IF (ln_cor) THEN
434            DO ji = 1, jnumcorio
435               WRITE(numout,'(1X,2A)') '             Coriolis input observation file name       coriofiles = ', &
436                  TRIM(coriofiles(ji))
437            END DO
438         ENDIF
439         IF (ln_profb) THEN
440            DO ji = 1, jnumprofb
441               IF (ln_profb_ena(ji)) THEN
442                  WRITE(numout,'(1X,2A)') '       Enact feedback input observation file name       profbfiles = ', &
443                     TRIM(profbfiles(ji))
444               ELSE
445                  WRITE(numout,'(1X,2A)') '             Feedback input observation file name       profbfiles = ', &
446                     TRIM(profbfiles(ji))
447               ENDIF
448               WRITE(numout,'(1X,2A)') '       Enact feedback input time setting switch    ln_profb_enatim = ', ln_profb_enatim(ji)
449            END DO
450         ENDIF
451         IF (ln_sladt) THEN
452            DO ji = 1, jnumslaact
453               WRITE(numout,'(1X,2A)') '             Active SLA input observation file name    slafilesact = ', &
454                  TRIM(slafilesact(ji))
455            END DO
456            DO ji = 1, jnumslapas
457               WRITE(numout,'(1X,2A)') '             Passive SLA input observation file name   slafilespas = ', &
458                  TRIM(slafilespas(ji))
459            END DO
460         ENDIF
461         IF (ln_slafb) THEN
462            DO ji = 1, jnumslafb
463               WRITE(numout,'(1X,2A)') '             Feedback SLA input observation file name   slafbfiles = ', &
464                  TRIM(slafbfiles(ji))
465            END DO
466         ENDIF
467         IF (ln_ghrsst) THEN
468            DO ji = 1, jnumsst
469               WRITE(numout,'(1X,2A)') '             GHRSST input observation file name           sstfiles = ', &
470                  TRIM(sstfiles(ji))
471            END DO
472         ENDIF
473         IF (ln_sstfb) THEN
474            DO ji = 1, jnumsstfb
475               WRITE(numout,'(1X,2A)') '             Feedback SST input observation file name   sstfbfiles = ', &
476                  TRIM(sstfbfiles(ji))
477            END DO
478         ENDIF
479         IF (ln_seaice) THEN
480            DO ji = 1, jnumseaice
481               WRITE(numout,'(1X,2A)') '             Sea Ice input observation file name       seaicefiles = ', &
482                  TRIM(seaicefiles(ji))
483            END DO
484         ENDIF
485         IF (ln_velavcur) THEN
486            DO ji = 1, jnumvelavcur
487               WRITE(numout,'(1X,2A)') '             Vel. cur. daily av. input file name     velavcurfiles = ', &
488                  TRIM(velavcurfiles(ji))
489            END DO
490         ENDIF
491         IF (ln_velhrcur) THEN
492            DO ji = 1, jnumvelhrcur
493               WRITE(numout,'(1X,2A)') '             Vel. cur. high freq. input file name    velhvcurfiles = ', &
494                  TRIM(velhrcurfiles(ji))
495            END DO
496         ENDIF
497         IF (ln_velavadcp) THEN
498            DO ji = 1, jnumvelavadcp
499               WRITE(numout,'(1X,2A)') '             Vel. ADCP daily av. input file name    velavadcpfiles = ', &
500                  TRIM(velavadcpfiles(ji))
501            END DO
502         ENDIF
503         IF (ln_velhradcp) THEN
504            DO ji = 1, jnumvelhradcp
505               WRITE(numout,'(1X,2A)') '             Vel. ADCP high freq. input file name   velhvadcpfiles = ', &
506                  TRIM(velhradcpfiles(ji))
507            END DO
508         ENDIF
509         IF (ln_velfb) THEN
510            DO ji = 1, jnumvelfb
511               IF (ln_velfb_av(ji)) THEN
512                  WRITE(numout,'(1X,2A)') '             Vel. feedback daily av. input file name    velfbfiles = ', &
513                     TRIM(velfbfiles(ji))
514               ELSE
515                  WRITE(numout,'(1X,2A)') '             Vel. feedback input observation file name  velfbfiles = ', &
516                     TRIM(velfbfiles(ji))
517               ENDIF
518            END DO
519         ENDIF
520         IF (ln_logchl) THEN
521            DO ji = 1, jnumlogchl
522               WRITE(numout,'(1X,2A)') '             logchl input observation file name        logchlfiles = ', &
523                  TRIM(logchlfiles(ji))
524            END DO
525         ENDIF
526         IF (ln_logchlfb) THEN
527            DO ji = 1, jnumlogchlfb
528               WRITE(numout,'(1X,2A)') '        Feedback logchl input observation file name  logchlfbfiles = ', &
529                  TRIM(logchlfbfiles(ji))
530            END DO
531         ENDIF
532         IF (ln_spm) THEN
533            DO ji = 1, jnumspm
534               WRITE(numout,'(1X,2A)') '             spm input observation file name              spmfiles = ', &
535                  TRIM(spmfiles(ji))
536            END DO
537         ENDIF
538         IF (ln_spmfb) THEN
539            DO ji = 1, jnumspmfb
540               WRITE(numout,'(1X,2A)') '             Feedback spm input observation file name   spmfbfiles = ', &
541                  TRIM(spmfbfiles(ji))
542            END DO
543         ENDIF
544         WRITE(numout,*) '             Initial date in window YYYYMMDD.HHMMSS        dobsini = ', dobsini
545         WRITE(numout,*) '             Final date in window YYYYMMDD.HHMMSS          dobsend = ', dobsend
546         WRITE(numout,*) '             Type of vertical interpolation method          n1dint = ', n1dint
547         WRITE(numout,*) '             Type of horizontal interpolation method        n2dint = ', n2dint
548         WRITE(numout,*) '             Rejection of observations near land swithch    ln_nea = ', ln_nea
549         WRITE(numout,*) '             MSSH correction scheme                         nmsshc = ', nmsshc
550         WRITE(numout,*) '             MDT  correction                               mdtcorr = ', mdtcorr
551         WRITE(numout,*) '             MDT cutoff for computed correction          mdtcutoff = ', mdtcutoff
552         WRITE(numout,*) '             Logical switch for alt bias                ln_altbias = ', ln_altbias
553         WRITE(numout,*) '             Logical switch for ignoring missing files   ln_ignmis = ', ln_ignmis
554         WRITE(numout,*) '             ENACT daily average types                             = ',endailyavtypes
555
556      ENDIF
557     
558      IF ( ln_vel3d .AND. ( .NOT. ln_grid_global ) ) THEN
559         CALL ctl_stop( 'Velocity data only works with ln_grid_global=.true.' )
560         RETURN
561      ENDIF
562
563      IF ( ln_grid_global ) THEN
564         CALL ctl_warn( 'ln_grid_global=T may cause memory issues when used with a large number of processors' )
565      ENDIF
566
567      CALL obs_typ_init
568     
569      IF ( ln_grid_global ) THEN     
570         CALL mppmap_init
571      ENDIF
572     
573      ! Parameter control
574#if defined key_diaobs
575      IF ( ( .NOT. ln_t3d ).AND.( .NOT. ln_s3d ).AND.( .NOT. ln_sla ).AND. &
576         & ( .NOT. ln_vel3d ).AND.                                         &
577         & ( .NOT. ln_ssh ).AND.( .NOT. ln_sst ).AND.( .NOT. ln_sss ).AND. &
578         & ( .NOT. ln_seaice ).AND.( .NOT. ln_vel3d ).AND.( .NOT. ln_logchl ).AND. &
579         & ( .NOT. ln_spm ) ) THEN
580         IF(lwp) WRITE(numout,cform_war)
581         IF(lwp) WRITE(numout,*) ' key_diaobs is activated but logical flags', &
582            &                    ' ln_t3d, ln_s3d, ln_sla, ln_ssh, ln_sst, ln_sss, ln_seaice, ln_vel3d,', &
583            &                    ' ln_logchl, ln_spm are all set to .FALSE.'
584         nwarn = nwarn + 1
585      ENDIF
586#endif
587
588      CALL obs_grid_setup( )
589      IF ( ( n1dint < 0 ).OR.( n1dint > 1 ) ) THEN
590         CALL ctl_stop(' Choice of vertical (1D) interpolation method', &
591            &                    ' is not available')
592      ENDIF
593      IF ( ( n2dint < 0 ).OR.( n2dint > 4 ) ) THEN
594         CALL ctl_stop(' Choice of horizontal (2D) interpolation method', &
595            &                    ' is not available')
596      ENDIF
597
598      !-----------------------------------------------------------------------
599      ! Depending on switches read the various observation types
600      !-----------------------------------------------------------------------
601      !  - Temperature/salinity profiles
602
603      IF ( ln_t3d .OR. ln_s3d ) THEN
604
605         ! Set the number of variables for profiles to 2 (T and S)
606         nprofvars = 2
607         ! Set the number of extra variables for profiles to 1 (insitu temp).
608         nprofextr = 1
609
610         ! Count how may insitu data sets we have and allocate data.
611         jprofset = 0
612         IF ( ln_ena ) jprofset = jprofset + 1
613         IF ( ln_cor ) jprofset = jprofset + 1
614         IF ( ln_profb ) jprofset = jprofset + jnumprofb
615         nprofsets = jprofset
616         IF ( nprofsets > 0 ) THEN
617            ALLOCATE(ld_enact(nprofsets))
618            ALLOCATE(profdata(nprofsets))
619            ALLOCATE(prodatqc(nprofsets))
620         ENDIF
621
622         jprofset = 0
623         
624         ! ENACT insitu data
625
626         IF ( ln_ena ) THEN
627
628            jprofset = jprofset + 1
629           
630            ld_enact(jprofset) = .TRUE.
631
632            CALL obs_rea_pro_dri( 1, profdata(jprofset),          &
633               &                  jnumenact, enactfiles(1:jnumenact), &
634               &                  nprofvars, nprofextr,        &
635               &                  nitend-nit000+2,             &
636               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
637               &                  ln_ignmis, ln_s_at_t, .TRUE., .FALSE., &
638               &                  kdailyavtypes = endailyavtypes )
639
640            DO jvar = 1, 2
641
642               CALL obs_prof_staend( profdata(jprofset), jvar )
643
644            END DO
645
646            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
647               &              ln_t3d, ln_s3d, ln_nea, &
648               &              kdailyavtypes=endailyavtypes )
649           
650         ENDIF
651
652         ! Coriolis insitu data
653
654         IF ( ln_cor ) THEN
655           
656            jprofset = jprofset + 1
657
658            ld_enact(jprofset) = .FALSE.
659
660            CALL obs_rea_pro_dri( 2, profdata(jprofset),          &
661               &                  jnumcorio, coriofiles(1:jnumcorio), &
662               &                  nprofvars, nprofextr,        &
663               &                  nitend-nit000+2,             &
664               &                  dobsini, dobsend, ln_t3d, ln_s3d, &
665               &                  ln_ignmis, ln_s_at_t, .FALSE., .FALSE. )
666
667            DO jvar = 1, 2
668
669               CALL obs_prof_staend( profdata(jprofset), jvar )
670
671            END DO
672
673            CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
674                 &            ln_t3d, ln_s3d, ln_nea )
675           
676         ENDIF
677 
678         ! Feedback insitu data
679
680         IF ( ln_profb ) THEN
681           
682            DO jset = 1, jnumprofb
683               
684               jprofset = jprofset + 1
685               ld_enact (jprofset) = ln_profb_ena(jset)
686
687               CALL obs_rea_pro_dri( 0, profdata(jprofset),          &
688                  &                  1, profbfiles(jset:jset), &
689                  &                  nprofvars, nprofextr,        &
690                  &                  nitend-nit000+2,             &
691                  &                  dobsini, dobsend, ln_t3d, ln_s3d, &
692                  &                  ln_ignmis, ln_s_at_t, &
693                  &                  ld_enact(jprofset).AND.&
694                  &                  ln_profb_enatim(jset), &
695                  &                  .FALSE., kdailyavtypes = endailyavtypes )
696               
697               DO jvar = 1, 2
698                 
699                  CALL obs_prof_staend( profdata(jprofset), jvar )
700                 
701               END DO
702               
703               IF ( ld_enact(jprofset) ) THEN
704                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
705                     &              ln_t3d, ln_s3d, ln_nea, &
706                     &              kdailyavtypes = endailyavtypes )
707               ELSE
708                  CALL obs_pre_pro( profdata(jprofset), prodatqc(jprofset),   &
709                     &              ln_t3d, ln_s3d, ln_nea )
710               ENDIF
711               
712            END DO
713
714         ENDIF
715
716      ENDIF
717
718      !  - Sea level anomalies
719      IF ( ln_sla ) THEN
720        ! Set the number of variables for sla to 1
721         nslavars = 1
722
723         ! Set the number of extra variables for sla to 2
724         nslaextr = 2
725         
726         ! Set the number of sla data sets to 2
727         nslasets = 0
728         IF ( ln_sladt ) THEN
729            nslasets = nslasets + 2
730         ENDIF
731         IF ( ln_slafb ) THEN
732            nslasets = nslasets + jnumslafb
733         ENDIF
734         
735         ALLOCATE(sladata(nslasets))
736         ALLOCATE(sladatqc(nslasets))
737         sladata(:)%nsurf=0
738         sladatqc(:)%nsurf=0
739
740         nslasets = 0
741
742         ! AVISO SLA data
743
744         IF ( ln_sladt ) THEN
745
746            ! Active SLA observations
747           
748            nslasets = nslasets + 1
749           
750            CALL obs_rea_sla( 1, sladata(nslasets), jnumslaact, &
751               &              slafilesact(1:jnumslaact), &
752               &              nslavars, nslaextr, nitend-nit000+2, &
753               &              dobsini, dobsend, ln_ignmis, .FALSE. )
754            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
755               &              ln_sla, ln_nea )
756           
757            ! Passive SLA observations
758           
759            nslasets = nslasets + 1
760           
761            CALL obs_rea_sla( 1, sladata(nslasets), jnumslapas, &
762               &              slafilespas(1:jnumslapas), &
763               &              nslavars, nslaextr, nitend-nit000+2, &
764               &              dobsini, dobsend, ln_ignmis, .FALSE. )
765           
766            CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
767               &              ln_sla, ln_nea )
768
769         ENDIF
770         
771         ! Feedback SLA data
772
773         IF ( ln_slafb ) THEN
774
775            DO jset = 1, jnumslafb
776           
777               nslasets = nslasets + 1
778           
779               CALL obs_rea_sla( 0, sladata(nslasets), 1, &
780                  &              slafbfiles(jset:jset), &
781                  &              nslavars, nslaextr, nitend-nit000+2, &
782                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
783               CALL obs_pre_sla( sladata(nslasets), sladatqc(nslasets), &
784                  &              ln_sla, ln_nea )
785
786            END DO               
787
788         ENDIF
789         
790         CALL obs_rea_mdt( nslasets, sladatqc, n2dint )
791           
792         ! read in altimeter bias
793         
794         IF ( ln_altbias ) THEN     
795            CALL obs_rea_altbias ( nslasets, sladatqc, n2dint, bias_file )
796         ENDIF
797     
798      ENDIF
799
800      !  - Sea surface height
801      IF ( ln_ssh ) THEN
802         IF(lwp) WRITE(numout,*) ' SSH currently not available'
803      ENDIF
804
805      !  - Sea surface temperature
806      IF ( ln_sst ) THEN
807
808         ! Set the number of variables for sst to 1
809         nsstvars = 1
810
811         ! Set the number of extra variables for sst to 0
812         nsstextr = 0
813
814         nsstsets = 0
815
816         IF (ln_reysst) nsstsets = nsstsets + 1
817         IF (ln_ghrsst) nsstsets = nsstsets + 1
818         IF ( ln_sstfb ) THEN
819            nsstsets = nsstsets + jnumsstfb
820         ENDIF
821
822         ALLOCATE(sstdata(nsstsets))
823         ALLOCATE(sstdatqc(nsstsets))
824         ALLOCATE(ld_sstnight(nsstsets))
825         sstdata(:)%nsurf=0
826         sstdatqc(:)%nsurf=0   
827         ld_sstnight(:)=.false.
828
829         nsstsets = 0
830
831         IF (ln_reysst) THEN
832
833            nsstsets = nsstsets + 1
834
835            ld_sstnight(nsstsets) = ln_sstnight
836
837            CALL obs_rea_sst_rey( reysstname, reysstfmt, sstdata(nsstsets), &
838               &                  nsstvars, nsstextr, &
839               &                  nitend-nit000+2, dobsini, dobsend )
840            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
841               &              ln_nea )
842
843        ENDIF
844       
845        IF (ln_ghrsst) THEN
846       
847            nsstsets = nsstsets + 1
848
849            ld_sstnight(nsstsets) = ln_sstnight
850         
851            CALL obs_rea_sst( 1, sstdata(nsstsets), jnumsst, &
852               &              sstfiles(1:jnumsst), &
853               &              nsstvars, nsstextr, nitend-nit000+2, &
854               &              dobsini, dobsend, ln_ignmis, .FALSE. )
855            CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), ln_sst, &
856               &              ln_nea )
857
858        ENDIF
859               
860         ! Feedback SST data
861
862         IF ( ln_sstfb ) THEN
863
864            DO jset = 1, jnumsstfb
865           
866               nsstsets = nsstsets + 1
867
868               ld_sstnight(nsstsets) = ln_sstnight
869           
870               CALL obs_rea_sst( 0, sstdata(nsstsets), 1, &
871                  &              sstfbfiles(jset:jset), &
872                  &              nsstvars, nsstextr, nitend-nit000+2, &
873                  &              dobsini, dobsend, ln_ignmis, .FALSE. )
874               CALL obs_pre_sst( sstdata(nsstsets), sstdatqc(nsstsets), &
875                  &              ln_sst, ln_nea )
876
877            END DO               
878
879         ENDIF
880         
881         !Read in bias field and correct SST.
882         IF ( ln_sstbias ) THEN
883            IF ( jnumsstbias == 0 ) CALL ctl_stop("ln_sstbias set,"// & 
884                                             "  but no bias"// & 
885                                             " files to read in")   
886            CALL obs_app_sstbias( nsstsets, sstdatqc, n2dint, & 
887                                  jnumsstbias, & 
888                                  sstbias_files(1:jnumsstbias) ) 
889         ENDIF
890
891      ENDIF
892
893      !  - Sea surface salinity
894      IF ( ln_sss ) THEN
895         IF(lwp) WRITE(numout,*) ' SSS currently not available'
896      ENDIF
897
898      !  - Sea Ice Concentration
899     
900      IF ( ln_seaice ) THEN
901
902         ! Set the number of variables for seaice to 1
903         nseaicevars = 1
904
905         ! Set the number of extra variables for seaice to 0
906         nseaiceextr = 0
907         
908         ! Set the number of data sets to 1
909         nseaicesets = 1
910
911         ALLOCATE(seaicedata(nseaicesets))
912         ALLOCATE(seaicedatqc(nseaicesets))
913         seaicedata(:)%nsurf=0
914         seaicedatqc(:)%nsurf=0
915
916         CALL obs_rea_seaice( 1, seaicedata(nseaicesets), jnumseaice, &
917            &                 seaicefiles(1:jnumseaice), &
918            &                 nseaicevars, nseaiceextr, nitend-nit000+2, &
919            &                 dobsini, dobsend, ln_ignmis, .FALSE. )
920
921         CALL obs_pre_seaice( seaicedata(nseaicesets), seaicedatqc(nseaicesets), &
922            &                 ln_seaice, ln_nea )
923 
924      ENDIF
925
926      IF (ln_vel3d) THEN
927
928         ! Set the number of variables for profiles to 2 (U and V)
929         nvelovars = 2
930
931         ! Set the number of extra variables for profiles to 2 to store
932         ! rotation parameters
933         nveloextr = 2
934
935         jveloset = 0
936         
937         IF ( ln_velavcur ) jveloset = jveloset + 1
938         IF ( ln_velhrcur ) jveloset = jveloset + 1
939         IF ( ln_velavadcp ) jveloset = jveloset + 1
940         IF ( ln_velhradcp ) jveloset = jveloset + 1
941         IF (ln_velfb) jveloset = jveloset + jnumvelfb
942
943         nvelosets = jveloset
944         IF ( nvelosets > 0 ) THEN
945            ALLOCATE( velodata(nvelosets) )
946            ALLOCATE( veldatqc(nvelosets) )
947            ALLOCATE( ld_velav(nvelosets) )
948         ENDIF
949         
950         jveloset = 0
951         
952         ! Daily averaged data
953
954         IF ( ln_velavcur ) THEN
955           
956            jveloset = jveloset + 1
957           
958            ld_velav(jveloset) = .TRUE.
959           
960            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavcur, &
961               &                  velavcurfiles(1:jnumvelavcur), &
962               &                  nvelovars, nveloextr, &
963               &                  nitend-nit000+2,              &
964               &                  dobsini, dobsend, ln_ignmis, &
965               &                  ld_velav(jveloset), &
966               &                  .FALSE. )
967           
968            DO jvar = 1, 2
969               CALL obs_prof_staend( velodata(jveloset), jvar )
970            END DO
971           
972            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
973               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
974           
975         ENDIF
976
977         ! High frequency data
978
979         IF ( ln_velhrcur ) THEN
980           
981            jveloset = jveloset + 1
982           
983            ld_velav(jveloset) = .FALSE.
984               
985            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhrcur, &
986               &                  velhrcurfiles(1:jnumvelhrcur), &
987               &                  nvelovars, nveloextr, &
988               &                  nitend-nit000+2,              &
989               &                  dobsini, dobsend, ln_ignmis, &
990               &                  ld_velav(jveloset), &
991               &                  .FALSE. )
992           
993            DO jvar = 1, 2
994               CALL obs_prof_staend( velodata(jveloset), jvar )
995            END DO
996           
997            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
998               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
999           
1000         ENDIF
1001
1002         ! Daily averaged data
1003
1004         IF ( ln_velavadcp ) THEN
1005           
1006            jveloset = jveloset + 1
1007           
1008            ld_velav(jveloset) = .TRUE.
1009           
1010            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelavadcp, &
1011               &                  velavadcpfiles(1:jnumvelavadcp), &
1012               &                  nvelovars, nveloextr, &
1013               &                  nitend-nit000+2,              &
1014               &                  dobsini, dobsend, ln_ignmis, &
1015               &                  ld_velav(jveloset), &
1016               &                  .FALSE. )
1017           
1018            DO jvar = 1, 2
1019               CALL obs_prof_staend( velodata(jveloset), jvar )
1020            END DO
1021           
1022            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1023               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1024           
1025         ENDIF
1026
1027         ! High frequency data
1028
1029         IF ( ln_velhradcp ) THEN
1030           
1031            jveloset = jveloset + 1
1032           
1033            ld_velav(jveloset) = .FALSE.
1034               
1035            CALL obs_rea_vel_dri( 1, velodata(jveloset), jnumvelhradcp, &
1036               &                  velhradcpfiles(1:jnumvelhradcp), &
1037               &                  nvelovars, nveloextr, &
1038               &                  nitend-nit000+2,              &
1039               &                  dobsini, dobsend, ln_ignmis, &
1040               &                  ld_velav(jveloset), &
1041               &                  .FALSE. )
1042           
1043            DO jvar = 1, 2
1044               CALL obs_prof_staend( velodata(jveloset), jvar )
1045            END DO
1046           
1047            CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1048               &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1049           
1050         ENDIF
1051
1052         IF ( ln_velfb ) THEN
1053
1054            DO jset = 1, jnumvelfb
1055           
1056               jveloset = jveloset + 1
1057
1058               ld_velav(jveloset) = ln_velfb_av(jset)
1059               
1060               CALL obs_rea_vel_dri( 0, velodata(jveloset), 1, &
1061                  &                  velfbfiles(jset:jset), &
1062                  &                  nvelovars, nveloextr, &
1063                  &                  nitend-nit000+2,              &
1064                  &                  dobsini, dobsend, ln_ignmis, &
1065                  &                  ld_velav(jveloset), &
1066                  &                  .FALSE. )
1067               
1068               DO jvar = 1, 2
1069                  CALL obs_prof_staend( velodata(jveloset), jvar )
1070               END DO
1071               
1072               CALL obs_pre_vel( velodata(jveloset), veldatqc(jveloset), &
1073                  &              ln_vel3d, ln_nea, ld_velav(jveloset) )
1074
1075
1076            END DO
1077           
1078         ENDIF
1079
1080      ENDIF
1081
1082      !  - log10(chlorophyll)
1083     
1084      IF ( ln_logchl ) THEN
1085
1086         ! Set the number of variables for logchl to 1
1087         nlogchlvars = 1
1088
1089         ! Set the number of extra variables for logchl to 0
1090         nlogchlextr = 0
1091         
1092         IF ( ln_logchlfb ) THEN
1093            nlogchlsets = jnumlogchlfb
1094         ELSE
1095            nlogchlsets = 1
1096         ENDIF
1097
1098         ALLOCATE(logchldata(nlogchlsets))
1099         ALLOCATE(logchldatqc(nlogchlsets))
1100         logchldata(:)%nsurf=0
1101         logchldatqc(:)%nsurf=0
1102
1103         nlogchlsets = 0
1104
1105         IF ( ln_logchlfb ) THEN             ! Feedback file format
1106
1107            DO jset = 1, jnumlogchlfb
1108           
1109               nlogchlsets = nlogchlsets + 1
1110
1111               CALL obs_rea_logchl( 0, logchldata(nlogchlsets), 1, &
1112                  &                 logchlfbfiles(jset:jset), &
1113                  &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
1114                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1115
1116               CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1117                  &                 ln_logchl, ln_nea )
1118           
1119            ENDDO
1120
1121         ELSE                              ! Original file format
1122
1123            nlogchlsets = nlogchlsets + 1
1124
1125            CALL obs_rea_logchl( 1, logchldata(nlogchlsets), jnumlogchl, &
1126               &                 logchlfiles(1:jnumlogchl), &
1127               &                 nlogchlvars, nlogchlextr, nitend-nit000+2, &
1128               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1129
1130            CALL obs_pre_logchl( logchldata(nlogchlsets), logchldatqc(nlogchlsets), &
1131               &                 ln_logchl, ln_nea )
1132
1133         ENDIF
1134 
1135      ENDIF
1136
1137      !  - spm
1138     
1139      IF ( ln_spm ) THEN
1140
1141         ! Set the number of variables for spm to 1
1142         nspmvars = 1
1143
1144         ! Set the number of extra variables for spm to 0
1145         nspmextr = 0
1146         
1147         IF ( ln_spmfb ) THEN
1148            nspmsets = jnumspmfb
1149         ELSE
1150            nspmsets = 1
1151         ENDIF
1152
1153         ALLOCATE(spmdata(nspmsets))
1154         ALLOCATE(spmdatqc(nspmsets))
1155         spmdata(:)%nsurf=0
1156         spmdatqc(:)%nsurf=0
1157
1158         nspmsets = 0
1159
1160         IF ( ln_spmfb ) THEN             ! Feedback file format
1161
1162            DO jset = 1, jnumspmfb
1163           
1164               nspmsets = nspmsets + 1
1165
1166               CALL obs_rea_spm( 0, spmdata(nspmsets), 1, &
1167                  &                 spmfbfiles(jset:jset), &
1168                  &                 nspmvars, nspmextr, nitend-nit000+2, &
1169                  &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1170
1171               CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1172                  &                 ln_spm, ln_nea )
1173           
1174            ENDDO
1175
1176         ELSE                              ! Original file format
1177
1178            nspmsets = nspmsets + 1
1179
1180            CALL obs_rea_spm( 1, spmdata(nspmsets), jnumspm, &
1181               &                 spmfiles(1:jnumspm), &
1182               &                 nspmvars, nspmextr, nitend-nit000+2, &
1183               &                 dobsini, dobsend, ln_ignmis, .FALSE. )
1184
1185            CALL obs_pre_spm( spmdata(nspmsets), spmdatqc(nspmsets), &
1186               &                 ln_spm, ln_nea )
1187
1188         ENDIF
1189 
1190      ENDIF
1191     
1192   END SUBROUTINE dia_obs_init
1193
1194   SUBROUTINE dia_obs( kstp )
1195      !!----------------------------------------------------------------------
1196      !!                    ***  ROUTINE dia_obs  ***
1197      !!         
1198      !! ** Purpose : Call the observation operators on each time step
1199      !!
1200      !! ** Method  : Call the observation operators on each time step to
1201      !!              compute the model equivalent of the following date:
1202      !!               - T profiles
1203      !!               - S profiles
1204      !!               - Sea surface height (referenced to a mean)
1205      !!               - Sea surface temperature
1206      !!               - Sea surface salinity
1207      !!               - Velocity component (U,V) profiles
1208      !!               - Sea surface log10(chlorophyll)
1209      !!               - Sea surface spm
1210      !!
1211      !! ** Action  :
1212      !!
1213      !! History :
1214      !!        !  06-03  (K. Mogensen) Original code
1215      !!        !  06-05  (K. Mogensen) Reformatted
1216      !!        !  06-10  (A. Weaver) Cleaning
1217      !!        !  07-03  (K. Mogensen) General handling of profiles
1218      !!        !  07-04  (G. Smith) Generalized surface operators
1219      !!        !  08-10  (M. Valdivieso) obs operator for velocity profiles
1220      !!        !  14-08  (J. While) observation operator for profiles in
1221      !!                             generalised vertical coordinates
1222      !!----------------------------------------------------------------------
1223      !! * Modules used
1224      USE dom_oce, ONLY : &             ! Ocean space and time domain variables
1225         & rdt,           &                       
1226         & gdept_1d,       &             
1227#if defined key_vvl 
1228         & gdept_n,       &
1229#else 
1230         & gdept_1d,      &
1231#endif                                       
1232         & tmask, umask, vmask                           
1233      USE phycst, ONLY : &              ! Physical constants
1234         & rday                         
1235      USE oce, ONLY : &                 ! Ocean dynamics and tracers variables
1236         & tsn,  &             
1237         & un, vn,  &
1238         & sshn
1239#if defined  key_lim3
1240      USE ice, ONLY : &                     ! LIM Ice model variables
1241         & frld
1242#endif
1243#if defined key_lim2
1244      USE ice_2, ONLY : &                     ! LIM Ice model variables
1245         & frld
1246#endif
1247#if defined key_hadocc
1248      USE trc, ONLY :  &                ! HadOCC chlorophyll
1249         & HADOCC_CHL, &
1250         & HADOCC_FILL_FLT
1251#elif defined key_medusa && defined key_foam_medusa
1252      USE trc, ONLY :  &                ! MEDUSA chlorophyll
1253         & MEDUSA_CHL, &
1254         & MEDUSA_FILL_FLT
1255#elif defined key_fabm
1256      !USE ???                           ! ERSEM chlorophyll
1257#endif
1258#if defined key_spm
1259      USE par_spm, ONLY: &              ! ERSEM/SPM sediments
1260         & jp_spm
1261#endif
1262      IMPLICIT NONE
1263
1264      !! * Arguments
1265      INTEGER, INTENT(IN) :: kstp                         ! Current timestep
1266      !! * Local declarations
1267      INTEGER :: idaystp                ! Number of timesteps per day
1268      INTEGER :: jprofset               ! Profile data set loop variable
1269      INTEGER :: jslaset                ! SLA data set loop variable
1270      INTEGER :: jsstset                ! SST data set loop variable
1271      INTEGER :: jseaiceset             ! sea ice data set loop variable
1272      INTEGER :: jveloset               ! velocity profile data loop variable
1273      INTEGER :: jlogchlset             ! logchl data set loop variable
1274      INTEGER :: jspmset                ! spm data set loop variable
1275      INTEGER :: jvar                   ! Variable number   
1276#if ! defined key_lim2 && ! defined key_lim3
1277      REAL(wp), POINTER, DIMENSION(:,:) :: frld   
1278#endif
1279      REAL(wp) :: tiny                  ! small number
1280      REAL(wp), DIMENSION(jpi,jpj) :: &
1281         logchl                         ! array for log chlorophyll
1282      REAL(wp), DIMENSION(jpi,jpj) :: &
1283         maskchl                        ! array for special chlorophyll mask
1284      REAL(wp), DIMENSION(jpi,jpj) :: &
1285         spm                            ! array for spm
1286      INTEGER :: jn                     ! loop index
1287      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1288 
1289#if ! defined key_lim2 && ! defined key_lim3
1290      CALL wrk_alloc(jpi,jpj,frld) 
1291#endif
1292
1293      IF(lwp) THEN
1294         WRITE(numout,*)
1295         WRITE(numout,*) 'dia_obs : Call the observation operators', kstp
1296         WRITE(numout,*) '~~~~~~~'
1297      ENDIF
1298
1299      idaystp = NINT( rday / rdt )
1300
1301      !-----------------------------------------------------------------------
1302      ! No LIM => frld == 0.0_wp
1303      !-----------------------------------------------------------------------
1304#if ! defined key_lim2 && ! defined key_lim3
1305      frld(:,:) = 0.0_wp
1306#endif
1307      !-----------------------------------------------------------------------
1308      ! Depending on switches call various observation operators
1309      !-----------------------------------------------------------------------
1310
1311      !  - Temperature/salinity profiles
1312      IF ( ln_t3d .OR. ln_s3d ) THEN
1313         DO jprofset = 1, nprofsets
1314            IF( (.NOT. lk_vvl) .AND. (ln_zco .OR. ln_zps) ) THEN
1315               IF(lwp) THEN
1316                  WRITE(numout,*) 'dia_obs : calling obs_pro_opt'
1317               ENDIF
1318               IF ( ld_enact(jprofset) ) THEN
1319                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1320                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1321                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1322                     &              gdept_1d, tmask, n1dint, n2dint,        & 
1323                     &              kdailyavtypes = endailyavtypes ) 
1324               ELSE
1325                  CALL obs_pro_opt( prodatqc(jprofset),                     & 
1326                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1327                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1328                     &              gdept_1d, tmask, n1dint, n2dint               ) 
1329               ENDIF
1330            ELSE
1331               IF(lwp) THEN
1332                  WRITE(numout,*) 'dia_obs : calling obs_pro_sco_opt'
1333               ENDIF
1334               IF ( ld_enact(jprofset) ) THEN
1335                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1336                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1337                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1338                     &              fsdept(:,:,:), fsdepw(:,:,:),           & 
1339                     &              tmask, n1dint, n2dint,                  & 
1340                     &              kdailyavtypes = endailyavtypes ) 
1341               ELSE
1342                  CALL obs_pro_sco_opt( prodatqc(jprofset),                 & 
1343                     &              kstp, jpi, jpj, jpk, nit000, idaystp,   & 
1344                     &              tsn(:,:,:,jp_tem), tsn(:,:,:,jp_sal),   & 
1345                     &              fsdept(:,:,:), fsdepw(:,:,:),           &
1346                     &              tmask, n1dint, n2dint ) 
1347               ENDIF
1348            ENDIF
1349         END DO
1350      ENDIF
1351
1352      !  - Sea surface anomaly
1353      IF ( ln_sla ) THEN
1354         DO jslaset = 1, nslasets
1355            CALL obs_sla_opt( sladatqc(jslaset),            &
1356               &              kstp, jpi, jpj, nit000, sshn, &
1357               &              tmask(:,:,1), n2dint )
1358         END DO         
1359      ENDIF
1360
1361      !  - Sea surface temperature
1362      IF ( ln_sst ) THEN
1363         DO jsstset = 1, nsstsets
1364            CALL obs_sst_opt( sstdatqc(jsstset),                &
1365               &              kstp, jpi, jpj, nit000, idaystp,  &
1366               &              tsn(:,:,1,jp_tem), tmask(:,:,1),  &
1367               &              n2dint, ld_sstnight(jsstset) )
1368         END DO
1369      ENDIF
1370
1371      !  - Sea surface salinity
1372      IF ( ln_sss ) THEN
1373         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1374      ENDIF
1375
1376#if defined key_lim2 || defined key_lim3
1377      IF ( ln_seaice ) THEN
1378         DO jseaiceset = 1, nseaicesets
1379            CALL obs_seaice_opt( seaicedatqc(jseaiceset),      &
1380               &              kstp, jpi, jpj, nit000, 1.-frld, &
1381               &              tmask(:,:,1), n2dint )
1382         END DO
1383      ENDIF     
1384#endif
1385
1386      !  - Velocity profiles
1387      IF ( ln_vel3d ) THEN
1388         DO jveloset = 1, nvelosets
1389           ! zonal component of velocity
1390           CALL obs_vel_opt( veldatqc(jveloset), kstp, jpi, jpj, jpk, &
1391              &              nit000, idaystp, un, vn, gdept_1d, umask, vmask, &
1392                             n1dint, n2dint, ld_velav(jveloset) )
1393         END DO
1394      ENDIF
1395
1396      IF ( ln_logchl ) THEN
1397
1398#if defined key_hadocc
1399         logchl(:,:)  = HADOCC_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
1400#elif defined key_medusa && defined key_foam_medusa
1401         logchl(:,:)  = MEDUSA_CHL(:,:,1)    ! (not log) chlorophyll from HadOCC
1402#elif defined key_fabm
1403         !logchl(:,:)  =  ???                 ! (not log) chlorophyll from ERSEM
1404         CALL ctl_stop( ' Trying to run logchl observation operator', &
1405            &           ' but not properly implemented for FABM-ERSEM yet' )
1406#else
1407         CALL ctl_stop( ' Trying to run logchl observation operator', &
1408            &           ' but no biogeochemical model appears to have been defined' )
1409#endif
1410
1411         maskchl(:,:) = tmask(:,:,1)         ! create a special mask to exclude certain things
1412
1413         ! Take the log10 where we can, otherwise exclude
1414         tiny = 1.0e-20
1415         WHERE(logchl(:,:) > tiny .AND. logchl(:,:) /= obfillflt )
1416            logchl(:,:)  = LOG10(logchl(:,:))
1417         ELSEWHERE
1418            logchl(:,:)  = obfillflt
1419            maskchl(:,:) = 0
1420         END WHERE
1421
1422         DO jlogchlset = 1, nlogchlsets
1423             CALL obs_logchl_opt( logchldatqc(jlogchlset),             &
1424               &                  kstp, jpi, jpj, nit000, logchl(:,:), &
1425               &                  maskchl(:,:), n2dint )
1426         END DO         
1427      ENDIF
1428
1429      IF ( ln_spm ) THEN
1430#if defined key_spm
1431         spm(:,:) = 0.0
1432         DO jn = 1, jp_spm
1433            spm(:,:) = spm(:,:) + trn(:,:,1,jn)   ! sum SPM sizes
1434         END DO
1435#else
1436         CALL ctl_stop( ' Trying to run spm observation operator', &
1437            &           ' but no spm model appears to have been defined' )
1438#endif
1439
1440         DO jspmset = 1, nspmsets
1441             CALL obs_spm_opt( spmdatqc(jspmset),                &
1442               &               kstp, jpi, jpj, nit000, spm(:,:), &
1443               &               tmask(:,:,1), n2dint )
1444         END DO         
1445      ENDIF
1446
1447#if ! defined key_lim2 && ! defined key_lim3
1448      CALL wrk_dealloc(jpi,jpj,frld) 
1449#endif
1450
1451   END SUBROUTINE dia_obs
1452 
1453   SUBROUTINE dia_obs_wri 
1454      !!----------------------------------------------------------------------
1455      !!                    ***  ROUTINE dia_obs_wri  ***
1456      !!         
1457      !! ** Purpose : Call observation diagnostic output routines
1458      !!
1459      !! ** Method  : Call observation diagnostic output routines
1460      !!
1461      !! ** Action  :
1462      !!
1463      !! History :
1464      !!        !  06-03  (K. Mogensen) Original code
1465      !!        !  06-05  (K. Mogensen) Reformatted
1466      !!        !  06-10  (A. Weaver) Cleaning
1467      !!        !  07-03  (K. Mogensen) General handling of profiles
1468      !!        !  08-09  (M. Valdivieso) Velocity component (U,V) profiles
1469      !!----------------------------------------------------------------------
1470      IMPLICIT NONE
1471
1472      !! * Local declarations
1473
1474      INTEGER :: jprofset                 ! Profile data set loop variable
1475      INTEGER :: jveloset                 ! Velocity data set loop variable
1476      INTEGER :: jslaset                  ! SLA data set loop variable
1477      INTEGER :: jsstset                  ! SST data set loop variable
1478      INTEGER :: jseaiceset               ! Sea Ice data set loop variable
1479      INTEGER :: jlogchlset               ! logchl data set loop variable
1480      INTEGER :: jspmset                  ! spm data set loop variable
1481      INTEGER :: jset
1482      INTEGER :: jfbini
1483      CHARACTER(LEN=20) :: datestr=" ",timestr=" "
1484      CHARACTER(LEN=10) :: cdtmp
1485      !-----------------------------------------------------------------------
1486      ! Depending on switches call various observation output routines
1487      !-----------------------------------------------------------------------
1488
1489      !  - Temperature/salinity profiles
1490
1491      IF( ln_t3d .OR. ln_s3d ) THEN
1492
1493         ! Copy data from prodatqc to profdata structures
1494         DO jprofset = 1, nprofsets
1495
1496            CALL obs_prof_decompress( prodatqc(jprofset), &
1497                 &                    profdata(jprofset), .TRUE., numout )
1498
1499         END DO
1500
1501         ! Write the profiles.
1502
1503         jprofset = 0
1504
1505         ! ENACT insitu data
1506
1507         IF ( ln_ena ) THEN
1508           
1509            jprofset = jprofset + 1
1510
1511            CALL obs_wri_p3d( 'enact', profdata(jprofset) )
1512
1513         ENDIF
1514
1515         ! Coriolis insitu data
1516
1517         IF ( ln_cor ) THEN
1518           
1519            jprofset = jprofset + 1
1520
1521            CALL obs_wri_p3d( 'corio', profdata(jprofset) )
1522           
1523         ENDIF
1524         
1525         ! Feedback insitu data
1526
1527         IF ( ln_profb ) THEN
1528
1529            jfbini = jprofset + 1
1530
1531            DO jprofset = jfbini, nprofsets
1532               
1533               jset = jprofset - jfbini + 1
1534               WRITE(cdtmp,'(A,I2.2)')'profb_',jset
1535               CALL obs_wri_p3d( cdtmp, profdata(jprofset) )
1536
1537            END DO
1538
1539         ENDIF
1540
1541      ENDIF
1542
1543      !  - Sea surface anomaly
1544      IF ( ln_sla ) THEN
1545
1546         ! Copy data from sladatqc to sladata structures
1547         DO jslaset = 1, nslasets
1548
1549              CALL obs_surf_decompress( sladatqc(jslaset), &
1550                 &                    sladata(jslaset), .TRUE., numout )
1551
1552         END DO
1553
1554         jslaset = 0 
1555
1556         ! Write the AVISO SLA data
1557
1558         IF ( ln_sladt ) THEN
1559           
1560            jslaset = 1
1561            CALL obs_wri_sla( 'aviso_act', sladata(jslaset) )
1562            jslaset = 2
1563            CALL obs_wri_sla( 'aviso_pas', sladata(jslaset) )
1564
1565         ENDIF
1566
1567         IF ( ln_slafb ) THEN
1568           
1569            jfbini = jslaset + 1
1570
1571            DO jslaset = jfbini, nslasets
1572               
1573               jset = jslaset - jfbini + 1
1574               WRITE(cdtmp,'(A,I2.2)')'slafb_',jset
1575               CALL obs_wri_sla( cdtmp, sladata(jslaset) )
1576
1577            END DO
1578
1579         ENDIF
1580
1581      ENDIF
1582
1583      !  - Sea surface temperature
1584      IF ( ln_sst ) THEN
1585
1586         ! Copy data from sstdatqc to sstdata structures
1587         DO jsstset = 1, nsstsets
1588     
1589              CALL obs_surf_decompress( sstdatqc(jsstset), &
1590                 &                    sstdata(jsstset), .TRUE., numout )
1591
1592         END DO
1593
1594         jsstset = 0 
1595
1596         ! Write the AVISO SST data
1597
1598         IF ( ln_reysst ) THEN
1599           
1600            jsstset = jsstset + 1
1601            CALL obs_wri_sst( 'reynolds', sstdata(jsstset) )
1602
1603         ENDIF
1604
1605         IF ( ln_ghrsst ) THEN
1606           
1607            jsstset = jsstset + 1
1608            CALL obs_wri_sst( 'ghr', sstdata(jsstset) )
1609
1610         ENDIF
1611
1612         IF ( ln_sstfb ) THEN
1613           
1614            jfbini = jsstset + 1
1615
1616            DO jsstset = jfbini, nsstsets
1617               
1618               jset = jsstset - jfbini + 1
1619               WRITE(cdtmp,'(A,I2.2)')'sstfb_',jset
1620               CALL obs_wri_sst( cdtmp, sstdata(jsstset) )
1621
1622            END DO
1623
1624         ENDIF
1625
1626      ENDIF
1627
1628      !  - Sea surface salinity
1629      IF ( ln_sss ) THEN
1630         IF(lwp) WRITE(numout,*) ' SSS currently not available'
1631      ENDIF
1632
1633      !  - Sea Ice Concentration
1634      IF ( ln_seaice ) THEN
1635
1636         ! Copy data from seaicedatqc to seaicedata structures
1637         DO jseaiceset = 1, nseaicesets
1638
1639              CALL obs_surf_decompress( seaicedatqc(jseaiceset), &
1640                 &                    seaicedata(jseaiceset), .TRUE., numout )
1641
1642         END DO
1643
1644         ! Write the Sea Ice data
1645         DO jseaiceset = 1, nseaicesets
1646     
1647            WRITE(cdtmp,'(A,I2.2)')'seaicefb_',jseaiceset
1648            CALL obs_wri_seaice( cdtmp, seaicedata(jseaiceset) )
1649
1650         END DO
1651
1652      ENDIF
1653     
1654      ! Velocity data
1655      IF( ln_vel3d ) THEN
1656
1657         ! Copy data from veldatqc to velodata structures
1658         DO jveloset = 1, nvelosets
1659
1660            CALL obs_prof_decompress( veldatqc(jveloset), &
1661                 &                    velodata(jveloset), .TRUE., numout )
1662
1663         END DO
1664
1665         ! Write the profiles.
1666
1667         jveloset = 0
1668
1669         ! Daily averaged data
1670
1671         IF ( ln_velavcur ) THEN
1672           
1673            jveloset = jveloset + 1
1674
1675            CALL obs_wri_vel( 'velavcurr', velodata(jveloset), n2dint )
1676
1677         ENDIF
1678
1679         ! High frequency data
1680
1681         IF ( ln_velhrcur ) THEN
1682           
1683            jveloset = jveloset + 1
1684
1685            CALL obs_wri_vel( 'velhrcurr', velodata(jveloset), n2dint )
1686
1687         ENDIF
1688
1689         ! Daily averaged data
1690
1691         IF ( ln_velavadcp ) THEN
1692           
1693            jveloset = jveloset + 1
1694
1695            CALL obs_wri_vel( 'velavadcp', velodata(jveloset), n2dint )
1696
1697         ENDIF
1698
1699         ! High frequency data
1700
1701         IF ( ln_velhradcp ) THEN
1702           
1703            jveloset = jveloset + 1
1704           
1705            CALL obs_wri_vel( 'velhradcp', velodata(jveloset), n2dint )
1706               
1707         ENDIF
1708
1709         ! Feedback velocity data
1710
1711         IF ( ln_velfb ) THEN
1712
1713            jfbini = jveloset + 1
1714
1715            DO jveloset = jfbini, nvelosets
1716               
1717               jset = jveloset - jfbini + 1
1718               WRITE(cdtmp,'(A,I2.2)')'velfb_',jset
1719               CALL obs_wri_vel( cdtmp, velodata(jveloset), n2dint )
1720
1721            END DO
1722
1723         ENDIF
1724         
1725      ENDIF
1726
1727      !  - log10(chlorophyll)
1728      IF ( ln_logchl ) THEN
1729
1730         ! Copy data from logchldatqc to logchldata structures
1731         DO jlogchlset = 1, nlogchlsets
1732
1733            CALL obs_surf_decompress( logchldatqc(jlogchlset), &
1734                 &                    logchldata(jlogchlset), .TRUE., numout )
1735
1736         END DO
1737         
1738         ! Mark as bad observations with no valid model counterpart due to activities in dia_obs
1739         ! Seem to need to set to fill value rather than marking as bad to be effective, so do both
1740         DO jlogchlset = 1, nlogchlsets
1741            WHERE ( logchldata(jlogchlset)%rmod(:,1) == obfillflt )
1742               logchldata(jlogchlset)%nqc(:)    = 1
1743               logchldata(jlogchlset)%robs(:,1) = obfillflt
1744            END WHERE
1745         END DO
1746
1747         ! Write the logchl data
1748         DO jlogchlset = 1, nlogchlsets
1749     
1750            WRITE(cdtmp,'(A,I2.2)')'logchlfb_',jlogchlset
1751            CALL obs_wri_logchl( cdtmp, logchldata(jlogchlset) )
1752
1753         END DO
1754
1755      ENDIF
1756
1757      !  - spm
1758      IF ( ln_spm ) THEN
1759
1760         ! Copy data from spmdatqc to spmdata structures
1761         DO jspmset = 1, nspmsets
1762
1763            CALL obs_surf_decompress( spmdatqc(jspmset), &
1764                 &                    spmdata(jspmset), .TRUE., numout )
1765
1766         END DO
1767
1768         ! Write the spm data
1769         DO jspmset = 1, nspmsets
1770     
1771            WRITE(cdtmp,'(A,I2.2)')'spmfb_',jspmset
1772            CALL obs_wri_spm( cdtmp, spmdata(jspmset) )
1773
1774         END DO
1775
1776      ENDIF
1777
1778   END SUBROUTINE dia_obs_wri
1779
1780   SUBROUTINE dia_obs_dealloc
1781      IMPLICIT NONE
1782      !!----------------------------------------------------------------------
1783      !!                    *** ROUTINE dia_obs_dealloc ***
1784      !!
1785      !!  ** Purpose : To deallocate data to enable the obs_oper online loop.
1786      !!               Specifically: dia_obs_init --> dia_obs --> dia_obs_wri
1787      !!
1788      !!  ** Method : Clean up various arrays left behind by the obs_oper.
1789      !!
1790      !!  ** Action :
1791      !!
1792      !!----------------------------------------------------------------------
1793      !! obs_grid deallocation
1794      CALL obs_grid_deallocate
1795
1796      !! diaobs deallocation
1797      IF ( nprofsets > 0 ) THEN
1798          DEALLOCATE(ld_enact, &
1799                  &  profdata, &
1800                  &  prodatqc)
1801      END IF
1802      IF ( ln_sla ) THEN
1803          DEALLOCATE(sladata, &
1804                  &  sladatqc)
1805      END IF
1806      IF ( ln_seaice ) THEN
1807          DEALLOCATE(sladata, &
1808                  &  sladatqc)
1809      END IF
1810      IF ( ln_sst ) THEN
1811          DEALLOCATE(sstdata, &
1812                  &  sstdatqc)
1813      END IF
1814      IF ( ln_vel3d ) THEN
1815          DEALLOCATE(ld_velav, &
1816                  &  velodata, &
1817                  &  veldatqc)
1818      END IF
1819   END SUBROUTINE dia_obs_dealloc
1820
1821   SUBROUTINE ini_date( ddobsini )
1822      !!----------------------------------------------------------------------
1823      !!                    ***  ROUTINE ini_date  ***
1824      !!         
1825      !! ** Purpose : Get initial data in double precision YYYYMMDD.HHMMSS format
1826      !!
1827      !! ** Method  : Get initial data in double precision YYYYMMDD.HHMMSS format
1828      !!
1829      !! ** Action  : Get initial data in double precision YYYYMMDD.HHMMSS format
1830      !!
1831      !! History :
1832      !!        !  06-03  (K. Mogensen)  Original code
1833      !!        !  06-05  (K. Mogensen)  Reformatted
1834      !!        !  06-10  (A. Weaver) Cleaning
1835      !!        !  06-10  (G. Smith) Calculates initial date the same as method for final date
1836      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1837      !!----------------------------------------------------------------------
1838      USE phycst, ONLY : &            ! Physical constants
1839         & rday
1840!      USE daymod, ONLY : &            ! Time variables
1841!         & nmonth_len           
1842      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1843         & rdt
1844
1845      IMPLICIT NONE
1846
1847      !! * Arguments
1848      REAL(KIND=dp), INTENT(OUT) :: ddobsini                         ! Initial date in YYYYMMDD.HHMMSS
1849
1850      !! * Local declarations
1851      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1852      INTEGER :: imon
1853      INTEGER :: iday
1854      INTEGER :: ihou
1855      INTEGER :: imin
1856      INTEGER :: imday         ! Number of days in month.
1857      REAL(KIND=wp) :: zdayfrc ! Fraction of day
1858
1859      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1860
1861      !!----------------------------------------------------------------------
1862      !! Initial date initialization (year, month, day, hour, minute)
1863      !! (This assumes that the initial date is for 00z))
1864      !!----------------------------------------------------------------------
1865      iyea =   ndate0 / 10000
1866      imon = ( ndate0 - iyea * 10000 ) / 100
1867      iday =   ndate0 - iyea * 10000 - imon * 100
1868      ihou = 0
1869      imin = 0
1870
1871      !!----------------------------------------------------------------------
1872      !! Compute number of days + number of hours + min since initial time
1873      !!----------------------------------------------------------------------
1874      iday = iday + ( nit000 -1 ) * rdt / rday
1875      zdayfrc = ( nit000 -1 ) * rdt / rday
1876      zdayfrc = zdayfrc - aint(zdayfrc)
1877      ihou = int( zdayfrc * 24 )
1878      imin = int( (zdayfrc * 24 - ihou) * 60 )
1879
1880      !!-----------------------------------------------------------------------
1881      !! Convert number of days (iday) into a real date
1882      !!----------------------------------------------------------------------
1883
1884      CALL calc_month_len( iyea, imonth_len )
1885     
1886      DO WHILE ( iday > imonth_len(imon) )
1887         iday = iday - imonth_len(imon)
1888         imon = imon + 1 
1889         IF ( imon > 12 ) THEN
1890            imon = 1
1891            iyea = iyea + 1
1892            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1893         ENDIF
1894      END DO
1895
1896      !!----------------------------------------------------------------------
1897      !! Convert it into YYYYMMDD.HHMMSS format.
1898      !!----------------------------------------------------------------------
1899      ddobsini = iyea * 10000_dp + imon * 100_dp + &
1900         &       iday + ihou * 0.01_dp + imin * 0.0001_dp
1901
1902
1903   END SUBROUTINE ini_date
1904
1905   SUBROUTINE fin_date( ddobsfin )
1906      !!----------------------------------------------------------------------
1907      !!                    ***  ROUTINE fin_date  ***
1908      !!         
1909      !! ** Purpose : Get final data in double precision YYYYMMDD.HHMMSS format
1910      !!
1911      !! ** Method  : Get final data in double precision YYYYMMDD.HHMMSS format
1912      !!
1913      !! ** Action  : Get final data in double precision YYYYMMDD.HHMMSS format
1914      !!
1915      !! History :
1916      !!        !  06-03  (K. Mogensen)  Original code
1917      !!        !  06-05  (K. Mogensen)  Reformatted
1918      !!        !  06-10  (A. Weaver) Cleaning
1919      !!        !  10-05  (D. Lea) Update to month length calculation for NEMO vn3.2
1920      !!----------------------------------------------------------------------
1921      USE phycst, ONLY : &            ! Physical constants
1922         & rday
1923!      USE daymod, ONLY : &            ! Time variables
1924!         & nmonth_len               
1925      USE dom_oce, ONLY : &           ! Ocean space and time domain variables
1926         & rdt
1927
1928      IMPLICIT NONE
1929
1930      !! * Arguments
1931      REAL(KIND=dp), INTENT(OUT) :: ddobsfin                   ! Final date in YYYYMMDD.HHMMSS
1932
1933      !! * Local declarations
1934      INTEGER :: iyea        ! date - (year, month, day, hour, minute)
1935      INTEGER :: imon
1936      INTEGER :: iday
1937      INTEGER :: ihou
1938      INTEGER :: imin
1939      INTEGER :: imday         ! Number of days in month.
1940      REAL(KIND=wp) :: zdayfrc       ! Fraction of day
1941         
1942      INTEGER, DIMENSION(12) ::   imonth_len    !: length in days of the months of the current year
1943           
1944      !-----------------------------------------------------------------------
1945      ! Initial date initialization (year, month, day, hour, minute)
1946      ! (This assumes that the initial date is for 00z)
1947      !-----------------------------------------------------------------------
1948      iyea =   ndate0 / 10000
1949      imon = ( ndate0 - iyea * 10000 ) / 100
1950      iday =   ndate0 - iyea * 10000 - imon * 100
1951      ihou = 0
1952      imin = 0
1953     
1954      !-----------------------------------------------------------------------
1955      ! Compute number of days + number of hours + min since initial time
1956      !-----------------------------------------------------------------------
1957      iday    = iday +  nitend  * rdt / rday
1958      zdayfrc =  nitend  * rdt / rday
1959      zdayfrc = zdayfrc - AINT( zdayfrc )
1960      ihou    = INT( zdayfrc * 24 )
1961      imin    = INT( ( zdayfrc * 24 - ihou ) * 60 )
1962
1963      !-----------------------------------------------------------------------
1964      ! Convert number of days (iday) into a real date
1965      !----------------------------------------------------------------------
1966
1967      CALL calc_month_len( iyea, imonth_len )
1968     
1969      DO WHILE ( iday > imonth_len(imon) )
1970         iday = iday - imonth_len(imon)
1971         imon = imon + 1 
1972         IF ( imon > 12 ) THEN
1973            imon = 1
1974            iyea = iyea + 1
1975            CALL calc_month_len( iyea, imonth_len )  ! update month lengths
1976         ENDIF
1977      END DO
1978
1979      !-----------------------------------------------------------------------
1980      ! Convert it into YYYYMMDD.HHMMSS format
1981      !-----------------------------------------------------------------------
1982      ddobsfin = iyea * 10000_dp + imon * 100_dp    + iday &
1983         &     + ihou * 0.01_dp  + imin * 0.0001_dp
1984
1985    END SUBROUTINE fin_date
1986   
1987END MODULE diaobs
Note: See TracBrowser for help on using the repository browser.