source: branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 3586

Last change on this file since 3586 was 3586, checked in by cbricaud, 8 years ago

add modification from dev_r3342_MERCATOR7_SST in dev_MERCATOR_2012_rev3555

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