source: trunk/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 9 years ago

Merge of 3.4beta into the trunk

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