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

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

gross simplification of stand alone observation operator

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