New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diaobs.F90 in branches/UKMO/dev_r4650_general_vert_coord_obsoper_removetides/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_removetides/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 7375

Last change on this file since 7375 was 7375, checked in by mattmartin, 5 years ago

Added code to branch to allow sla model counterpart to be a time-average.

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