source: branches/UKMO/dev_r5518_seaiceobs/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6283

Last change on this file since 6283 was 6283, checked in by drew, 6 years ago

Commit changes Jennie made at VN3.5.

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