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.
obs_read_sla.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sla.F90 @ 4428

Last change on this file since 4428 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 20.6 KB
Line 
1MODULE obs_read_sla
2   !!======================================================================
3   !!                       ***  MODULE obs_read_sla  ***
4   !! Observation diagnostics: Read the along track SLA data from
5   !!                          the AVISO/CLS database or any SLA data
6   !!                          from feedback files
7   !!======================================================================
8
9   !!----------------------------------------------------------------------
10   !!   obs_rea_sla : Driver for reading SLA data from the AVISO/CLS/feedback
11   !!----------------------------------------------------------------------
12
13   !! * Modules used   
14   USE par_kind                 ! Precision variables
15   USE in_out_manager           ! I/O manager
16   USE dom_oce                  ! Ocean space and time domain variables
17   USE obs_mpp                  ! MPP support routines for observation diagnostics
18   USE julian                   ! Julian date routines
19   USE obs_utils                ! Observation operator utility functions
20   USE obs_grid                 ! Grid search
21   USE obs_sort                 ! Sorting observation arrays
22   USE obs_surf_def             ! Surface observation definitions
23   USE obs_types                ! Observation type definitions
24   USE obs_sla_io               ! I/O for sla files
25   USE netcdf                   ! NetCDF library
26
27   IMPLICIT NONE
28
29   !! * Routine accessibility
30   PRIVATE
31
32   PUBLIC obs_rea_sla  ! Read the SLA observations from the AVISO/SLA database
33
34   !! * Control permutation of array indices
35#  include "dom_oce_ftrans.h90"
36
37   !!----------------------------------------------------------------------
38   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
39   !! $Id$
40   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
41   !!----------------------------------------------------------------------
42
43CONTAINS
44
45   SUBROUTINE obs_rea_sla( kformat, &
46      &                    sladata, knumfiles, cfilenames, &
47      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
48      &                    ldignmis, ldmod, ldobstd )
49      !!---------------------------------------------------------------------
50      !!
51      !!                   *** ROUTINE obs_rea_sla ***
52      !!
53      !! ** Purpose : Read from file the SLA data
54      !!
55      !! ** Method  : Depending on kformat either AVISO or
56      !!              feedback data files are read
57      !!
58      !! ** Action  :
59      !!
60      !!
61      !! History : 
62      !!      ! :  2009-01 (K. Mogensen) Initial version based on old versions
63      !!----------------------------------------------------------------------
64      !! * Modules used
65
66      !! * Arguments
67      INTEGER :: kformat    ! Format of input data
68      !                     ! 0: Feedback
69      !                     ! 1: AVISO
70      TYPE(obs_surf), INTENT(INOUT) :: sladata    ! SLA data to be read
71      INTEGER, INTENT(IN) :: knumfiles      ! Number of files to read in
72      CHARACTER(LEN=128), INTENT(IN) ::  &
73         & cfilenames(knumfiles) ! File names to read in
74      INTEGER, INTENT(IN) :: kvars     ! Number of variables in sladata
75      INTEGER, INTENT(IN) :: kextr     ! Number of extra fields for each var in sladata
76      INTEGER, INTENT(IN) :: kstp      ! Ocean time-step index
77      LOGICAL, INTENT(IN) :: ldignmis    ! Ignore missing files
78      LOGICAL, INTENT(IN) :: ldmod       ! Initialize model from input data
79      LOGICAL, INTENT(INOUT), optional :: &
80         & ldobstd        ! Read observation standard deviation from fb. file
81      REAL(KIND=dp), INTENT(IN) :: ddobsini  ! Obs. ini time in YYYYMMDD.HHMMSS
82      REAL(KIND=dp), INTENT(IN) :: ddobsend  ! Obs. end time in YYYYMMDD.HHMMSS
83         
84      !! * Local declarations
85      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_sla'
86      INTEGER :: ji
87      INTEGER :: jj
88      INTEGER :: jk
89      INTEGER :: iflag
90      INTEGER :: inobf
91      INTEGER :: i_file_id
92      INTEGER :: inowin
93      INTEGER :: iyea
94      INTEGER :: imon
95      INTEGER :: iday
96      INTEGER :: ihou
97      INTEGER :: imin
98      INTEGER :: isec
99      INTEGER, DIMENSION(knumfiles) :: &
100         & irefdate
101      INTEGER :: iobsmpp
102      INTEGER, DIMENSION(imaxmissions+1) :: &
103         & ityp, &
104         & itypmpp
105      INTEGER, DIMENSION(:), ALLOCATABLE :: &
106         & iobsi,    &
107         & iobsj,    &
108         & iproc,    &
109         & iindx,    &
110         & ifileidx, &
111         & islaidx
112      INTEGER :: itype 
113      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
114         & zphi, &
115         & zlam
116      REAL(dp), DIMENSION(:), ALLOCATABLE :: &
117         & zdat
118      LOGICAL :: llvalprof
119      LOGICAL :: llobstd
120      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
121         & inpfiles
122      REAL(dp), DIMENSION(knumfiles) :: &
123         & djulini, &
124         & djulend
125      INTEGER, DIMENSION(knumfiles) :: &
126         & iobspos, &
127         & iobcpos
128      INTEGER :: iobs
129      INTEGER :: iobstot
130      CHARACTER(len=8) :: cl_refdate
131   
132      ! Local initialization
133      iobs = 0
134      IF ( PRESENT(ldobstd) ) THEN
135         IF (.NOT.ldmod) THEN
136            llobstd = .false.
137         ELSE
138            llobstd = ldobstd
139         ENDIF
140      ELSE
141         llobstd = .FALSE.
142      ENDIF
143 
144      !-----------------------------------------------------------------------
145      ! Check that the model part is just with feedback data files
146      !-----------------------------------------------------------------------
147      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
148         CALL ctl_stop( 'Model can only be read from feedback data' )
149         RETURN
150      ENDIF
151
152      !-----------------------------------------------------------------------
153      ! Check that the prescribed obs err is just with feedback data files
154      !-----------------------------------------------------------------------
155      IF ( llobstd .AND. ( kformat /= 0 ) ) THEN
156         CALL ctl_stop( 'Observation error can only be read from feedback files' )
157         RETURN
158      ENDIF
159
160      !-----------------------------------------------------------------------
161      ! Count the number of files needed and allocate the obfbdata type
162      !-----------------------------------------------------------------------
163     
164      inobf = knumfiles
165     
166      ALLOCATE( inpfiles(inobf) )
167
168      sla_files : DO jj = 1, inobf
169         
170         !---------------------------------------------------------------------
171         ! Prints
172         !---------------------------------------------------------------------
173         IF(lwp) THEN
174            WRITE(numout,*)
175            WRITE(numout,*) ' obs_rea_sla : Reading from file = ', &
176               & TRIM( TRIM( cfilenames(jj) ) )
177            WRITE(numout,*) ' ~~~~~~~~~~~'
178            WRITE(numout,*)
179         ENDIF
180
181         !---------------------------------------------------------------------
182         !  Initialization: Open file and get dimensions only
183         !---------------------------------------------------------------------
184         
185         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
186            &                      i_file_id )
187         
188         IF ( iflag /= nf90_noerr ) THEN
189
190            IF ( ldignmis ) THEN
191               inpfiles(jj)%nobs = 0
192               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
193                  &           ' not found' )
194            ELSE
195               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
196                  &           ' not found' )
197            ENDIF
198
199         ELSE 
200           
201            !------------------------------------------------------------------
202            !  Close the file since it is opened in read_proffile
203            !------------------------------------------------------------------
204           
205            iflag = nf90_close( i_file_id )
206
207            !------------------------------------------------------------------
208            !  Read the profile file into inpfiles
209            !------------------------------------------------------------------
210            IF ( kformat == 0 ) THEN
211               CALL init_obfbdata( inpfiles(jj) )
212               IF(lwp) THEN
213                  WRITE(numout,*)
214                  WRITE(numout,*)'Reading from feedback file :', &
215                     &           TRIM( cfilenames(jj) )
216               ENDIF
217               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
218                  &                ldgrid = .TRUE. )
219               IF ( inpfiles(jj)%nvar < 1 ) THEN
220                  CALL ctl_stop( 'Feedback format error' )
221                  RETURN
222               ENDIF
223               IF ( TRIM(inpfiles(jj)%cname(1)) /= 'SLA' ) THEN
224                  CALL ctl_stop( 'Feedback format error' )
225                  RETURN
226               ENDIF
227               IF ( ldmod .AND. ( ( inpfiles(jj)%nadd == 0 ) .OR.&
228                  &               ( inpfiles(jj)%next == 0 ) ) ) THEN
229                  CALL ctl_stop( 'Model not in input data' )
230                  RETURN
231               ENDIF
232            ELSEIF ( kformat == 1) THEN
233               CALL read_avisofile( TRIM( cfilenames(jj) ), inpfiles(jj), &
234               &                    numout, lwp, .TRUE. )
235            ELSE
236               CALL ctl_stop( 'File format unknown' )
237            ENDIF
238
239            !------------------------------------------------------------------
240            !  Change longitude (-180,180)
241            !------------------------------------------------------------------
242
243            DO ji = 1, inpfiles(jj)%nobs 
244
245               IF ( inpfiles(jj)%plam(ji) < -180. ) &
246                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
247
248               IF ( inpfiles(jj)%plam(ji) >  180. ) &
249                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
250
251            END DO
252
253            !------------------------------------------------------------------
254            !  Calculate the date  (change eventually)
255            !------------------------------------------------------------------
256            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
257            READ(cl_refdate,'(I8)') irefdate(jj)
258           
259            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
260            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
261               &           krefdate = irefdate(jj) )
262            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
263            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
264               &           krefdate = irefdate(jj) )
265            IF ( inpfiles(jj)%nobs > 0 ) THEN
266               inpfiles(jj)%iproc = -1
267               inpfiles(jj)%iobsi = -1
268               inpfiles(jj)%iobsj = -1
269            ENDIF
270            inowin = 0
271            DO ji = 1, inpfiles(jj)%nobs
272               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
273               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
274               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
275                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
276                  inowin = inowin + 1
277               ENDIF
278            END DO
279            ALLOCATE( zlam(inowin)  )
280            ALLOCATE( zphi(inowin)  )
281            ALLOCATE( iobsi(inowin) )
282            ALLOCATE( iobsj(inowin) )
283            ALLOCATE( iproc(inowin) )
284            inowin = 0
285            DO ji = 1, inpfiles(jj)%nobs
286               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
287               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
288               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
289                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
290                  inowin = inowin + 1
291                  zlam(inowin) = inpfiles(jj)%plam(ji)
292                  zphi(inowin) = inpfiles(jj)%pphi(ji)
293               ENDIF
294            END DO
295
296            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
297
298            inowin = 0
299            DO ji = 1, inpfiles(jj)%nobs
300               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
301               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
302               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
303                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
304                  inowin = inowin + 1
305                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
306                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
307                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
308               ENDIF
309            END DO
310            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
311
312            DO ji = 1, inpfiles(jj)%nobs
313               IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
314               IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
315               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
316                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
317                  IF ( nproc == 0 ) THEN
318                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
319                  ELSE
320                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
321                  ENDIF
322                  llvalprof = .FALSE.
323                  IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
324                     & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
325                     iobs = iobs + 1
326                  ENDIF
327               ENDIF
328            END DO
329
330         ENDIF
331     
332      END DO sla_files
333
334      IF (llobstd) THEN
335         
336         DO jj = 1, inobf
337            iobspos(jj) = -1
338            iobcpos(jj) = -1
339            DO ji = 1,inpfiles(jj)%nadd
340               IF ( TRIM(inpfiles(jj)%caddname(ji)) == 'OSTD' ) THEN
341                  iobspos(jj)=ji
342               ENDIF
343               IF ( TRIM(inpfiles(jj)%caddname(ji)) == 'OCNT' ) THEN
344                  iobcpos(jj)=ji
345               ENDIF
346            END DO
347         END DO
348         llobstd = ( ( MINVAL(iobspos) > 0 ) .AND. ( MINVAL(iobcpos) > 0 ) )
349         IF (llobstd) THEN
350            IF (lwp) WRITE(numout,*)'SLA superobs information present.'
351         ELSE
352            IF (lwp) WRITE(numout,*)'SLA superobs information not present.'
353         ENDIF
354
355      ENDIF
356     
357      !-----------------------------------------------------------------------
358      ! Get the time ordered indices of the input data
359      !-----------------------------------------------------------------------
360
361      !---------------------------------------------------------------------
362      !  Loop over input data files to count total number of profiles
363      !---------------------------------------------------------------------
364      iobstot = 0
365      DO jj = 1, inobf
366         DO ji = 1, inpfiles(jj)%nobs
367            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
368            IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
369            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
370               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
371               iobstot = iobstot + 1
372            ENDIF
373         END DO
374      END DO
375
376      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
377         &      islaidx(iobstot), zdat(iobstot) )
378      jk = 0
379      DO jj = 1, inobf
380         DO ji = 1, inpfiles(jj)%nobs
381            IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
382            IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
383            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
384               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
385               jk = jk + 1
386               ifileidx(jk) = jj
387               islaidx(jk) = ji
388               zdat(jk)     = inpfiles(jj)%ptim(ji)
389            ENDIF
390         END DO
391      END DO
392      CALL sort_dp_indx( iobstot, &
393         &               zdat,     &
394         &               iindx   )
395     
396      CALL obs_surf_alloc( sladata, iobs, kvars, kextr, kstp )
397     
398      ! * Read obs/positions, QC, all variable and assign to sladata
399 
400      iobs = 0
401
402      ityp   (:) = 0
403      itypmpp(:) = 0
404     
405      DO jk = 1, iobstot
406         
407         jj = ifileidx(iindx(jk))
408         ji = islaidx(iindx(jk))
409
410         IF ( inpfiles(jj)%ioqc(ji) > 2 ) CYCLE
411         IF ( inpfiles(jj)%ivqc(ji,1) > 2 ) CYCLE
412
413         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
414            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
415           
416            IF ( nproc == 0 ) THEN
417               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
418            ELSE
419               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
420            ENDIF
421           
422            ! Set observation information
423           
424            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
425               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
426
427               iobs = iobs + 1
428
429               CALL jul2greg( isec,                   &
430                  &           imin,                   &
431                  &           ihou,                   &
432                  &           iday,                   &
433                  &           imon,                   &
434                  &           iyea,                   &
435                  &           inpfiles(jj)%ptim(ji), &
436                  &           irefdate(jj) )
437
438
439               ! SLA time coordinates
440               sladata%nyea(iobs) = iyea
441               sladata%nmon(iobs) = imon
442               sladata%nday(iobs) = iday
443               sladata%nhou(iobs) = ihou
444               sladata%nmin(iobs) = imin
445               
446               ! SLA space coordinates
447               sladata%rlam(iobs) = inpfiles(jj)%plam(ji)
448               sladata%rphi(iobs) = inpfiles(jj)%pphi(ji)
449
450               ! Coordinate search parameters
451               sladata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
452               sladata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
453               
454               ! Instrument type
455               READ( inpfiles(jj)%cdtyp(ji), '(I4)' ) itype
456               sladata%ntyp(iobs) = itype
457               ityp(itype+1) = ityp(itype+1) + 1
458
459               ! Identifier
460               sladata%cwmo(iobs) = inpfiles(jj)%cdwmo(ji)
461
462               ! Bookkeeping data to match observations
463               sladata%nsidx(iobs) = iobs
464               sladata%nsfil(iobs) = iindx(jk)
465
466               ! QC flags
467               sladata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
468
469               ! Observed value
470               sladata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
471
472
473               ! Model and MDT is set to fbrmdi unless read from file
474               IF ( ldmod ) THEN
475                  sladata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
476                  sladata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
477                  sladata%rext(iobs,2) = inpfiles(jj)%pext(1,ji,1)
478                  IF (llobstd) THEN
479                     sladata%rext(iobs,3) = &
480                        & inpfiles(jj)%padd(1,ji,iobspos(jj),1)
481                     sladata%rext(iobs,4) = &
482                        & inpfiles(jj)%padd(1,ji,iobcpos(jj),1)
483                  ENDIF
484               ELSE
485                  sladata%rmod(iobs,1) = fbrmdi
486                  sladata%rext(iobs,:) = fbrmdi
487               ENDIF
488
489            ENDIF
490         ENDIF
491
492      END DO
493
494      !-----------------------------------------------------------------------
495      ! Sum up over processors
496      !-----------------------------------------------------------------------
497     
498      CALL obs_mpp_sum_integer( iobs, iobsmpp )
499      CALL obs_mpp_sum_integers( ityp, itypmpp, imaxmissions + 1 )
500     
501      !-----------------------------------------------------------------------
502      ! Output number of observations.
503      !-----------------------------------------------------------------------
504      IF (lwp) THEN
505
506         WRITE(numout,*)
507         WRITE(numout,'(1X,A)')'Altimeter satellites'
508         WRITE(numout,'(1X,A)')'--------------------'
509         DO jj = 1,8
510            IF ( itypmpp(jj) > 0 ) THEN
511               WRITE(numout,'(1X,A38,A2,I10)')calttyp(jj-1),'= ',itypmpp(jj)
512            ENDIF
513         END DO
514         WRITE(numout,'(1X,A50)')'--------------------------------------------------'
515         WRITE(numout,'(1X,A40,I10)')'Total                                 = ',iobsmpp
516         WRITE(numout,*)
517
518      ENDIF
519
520      !-----------------------------------------------------------------------
521      ! Deallocate temporary data
522      !-----------------------------------------------------------------------
523      DEALLOCATE( ifileidx, islaidx, zdat )
524
525      !-----------------------------------------------------------------------
526      ! Deallocate input data
527      !-----------------------------------------------------------------------
528      DO jj = 1, inobf
529         CALL dealloc_obfbdata( inpfiles(jj) )
530      END DO
531      DEALLOCATE( inpfiles )
532
533      !-----------------------------------------------------------------------
534      ! Reset ldobstd if the data is present
535      !-----------------------------------------------------------------------
536      IF ( PRESENT(ldobstd) ) THEN
537         ldobstd = llobstd
538      ENDIF
539
540   END SUBROUTINE obs_rea_sla
541
542END MODULE obs_read_sla
Note: See TracBrowser for help on using the repository browser.