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/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90 @ 6990

Last change on this file since 6990 was 6990, checked in by kingr, 6 years ago

Added code from nemo3.4 OBS branch to allow rejection of observations near open boundaries.

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