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_surf.F90 in branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_surf.F90 @ 5682

Last change on this file since 5682 was 5682, checked in by mattmartin, 9 years ago

OBS simplification changes committed to branch after running SETTE tests to make sure we get the same results as the trunk for ORCA2_LIM_OBS.
File size: 18.7 KB
Line 
1MODULE obs_read_surf
2   !!======================================================================
3   !!                       ***  MODULE obs_read_surf  ***
4   !! Observation diagnostics: Read the surface data from feedback files
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_rea_surf : Driver for reading surface data from feedback files
9   !!----------------------------------------------------------------------
10
11   !! * Modules used
12   USE par_kind                 ! Precision variables
13   USE in_out_manager           ! I/O manager
14   USE dom_oce                  ! Ocean space and time domain variables
15   USE obs_mpp                  ! MPP support routines for observation diagnostics
16   USE julian                   ! Julian date routines
17   USE obs_utils                ! Observation operator utility functions
18   USE obs_grid                 ! Grid search
19   USE obs_sort                 ! Sorting observation arrays
20   USE obs_surf_def             ! Surface observation definitions
21   USE obs_types                ! Observation type definitions
22   USE obs_fbm                  ! Feedback routines
23   USE netcdf                   ! NetCDF library
24
25   IMPLICIT NONE
26
27   !! * Routine accessibility
28   PRIVATE
29
30   PUBLIC obs_rea_surf      ! Read the surface observations from the point data
31
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
34   !! $Id: obs_read_surf.F90 4990 2014-12-15 16:42:49Z timgraham $
35   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37
38CONTAINS
39
40   SUBROUTINE obs_rea_surf( surfdata, knumfiles, cdfilenames, &
41      &                     kvars, kextr, kstp, ddobsini, ddobsend, &
42      &                     ldignmis, ldmod, ldnightav )
43      !!---------------------------------------------------------------------
44      !!
45      !!                   *** ROUTINE obs_rea_surf ***
46      !!
47      !! ** Purpose : Read from file the surface data
48      !!
49      !! ** Method  : Read in the data from feedback format files and
50      !!              put into the NEMO internal surface data structure
51      !!
52      !! ** Action  :
53      !!
54      !!
55      !! History : 
56      !!      ! :  2009-01 (K. Mogensen) Initial version based on old versions
57      !!      ! :  2015-02 (M. Martin)   Unify the different surface data type reading.
58      !!----------------------------------------------------------------------
59      !! * Modules used
60
61      !! * Arguments
62      TYPE(obs_surf), INTENT(INOUT) :: &
63         & surfdata                     ! Surface data to be read
64      INTEGER, INTENT(IN) :: knumfiles  ! Number of corio format files to read
65      CHARACTER(LEN=128), INTENT(IN) :: &
66         & cdfilenames(knumfiles)       ! File names to read in
67      INTEGER, INTENT(IN) :: kvars      ! Number of variables in surfdata
68      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var
69      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
70      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
71      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
72      LOGICAL, INTENT(IN) :: ldnightav  ! Observations represent a night-time average
73      REAL(dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
74      REAL(dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
75
76      !! * Local declarations
77      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_surf'
78      CHARACTER(len=8) :: clrefdate
79      CHARACTER(len=6), DIMENSION(:), ALLOCATABLE :: clvars
80      INTEGER :: ji
81      INTEGER :: jj
82      INTEGER :: jk
83      INTEGER :: iflag
84      INTEGER :: inobf
85      INTEGER :: i_file_id
86      INTEGER :: inowin
87      INTEGER :: iyea
88      INTEGER :: imon
89      INTEGER :: iday
90      INTEGER :: ihou
91      INTEGER :: imin
92      INTEGER :: isec
93      INTEGER :: itype
94      INTEGER :: iobsmpp
95      INTEGER :: iobs
96      INTEGER :: iobstot
97      INTEGER :: ios
98      INTEGER :: ioserrcount
99      INTEGER, PARAMETER :: jpsurfmaxtype = 1024
100      INTEGER, DIMENSION(knumfiles) :: irefdate
101      INTEGER, DIMENSION(jpsurfmaxtype+1) :: &
102         & ityp, &
103         & itypmpp
104      INTEGER, DIMENSION(:), ALLOCATABLE :: &
105         & iobsi,    &
106         & iobsj,    &
107         & iproc,    &
108         & iindx,    &
109         & ifileidx, &
110         & isurfidx
111      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
112         & zphi, &
113         & zlam
114      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
115         & zdat
116      REAL(wp), DIMENSION(knumfiles) :: &
117         & djulini, &
118         & djulend
119      LOGICAL :: llvalprof
120      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
121         & inpfiles
122
123      ! Local initialization
124      iobs = 0
125
126      !-----------------------------------------------------------------------
127      ! Count the number of files needed and allocate the obfbdata type
128      !-----------------------------------------------------------------------
129
130      inobf = knumfiles
131
132      ALLOCATE( inpfiles(inobf) )
133
134      surf_files : DO jj = 1, inobf
135
136         !---------------------------------------------------------------------
137         ! Prints
138         !---------------------------------------------------------------------
139         IF(lwp) THEN
140            WRITE(numout,*)
141            WRITE(numout,*) ' obs_rea_surf : Reading from file = ', &
142               & TRIM( TRIM( cdfilenames(jj) ) )
143            WRITE(numout,*) ' ~~~~~~~~~~~'
144            WRITE(numout,*)
145         ENDIF
146
147         !---------------------------------------------------------------------
148         !  Initialization: Open file and get dimensions only
149         !---------------------------------------------------------------------
150
151         iflag = nf90_open( TRIM( TRIM( cdfilenames(jj) ) ), nf90_nowrite, &
152            &                      i_file_id )
153
154         IF ( iflag /= nf90_noerr ) THEN
155
156            IF ( ldignmis ) THEN
157               inpfiles(jj)%nobs = 0
158               CALL ctl_warn( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
159                  &           ' not found' )
160            ELSE
161               CALL ctl_stop( 'File ' // TRIM( TRIM( cdfilenames(jj) ) ) // &
162                  &           ' not found' )
163            ENDIF
164
165         ELSE 
166
167            !------------------------------------------------------------------
168            !  Close the file since it is opened in read_obfbdata
169            !------------------------------------------------------------------
170
171            iflag = nf90_close( i_file_id )
172
173            !------------------------------------------------------------------
174            !  Read the profile file into inpfiles
175            !------------------------------------------------------------------
176            CALL init_obfbdata( inpfiles(jj) )
177            CALL read_obfbdata( TRIM( cdfilenames(jj) ), inpfiles(jj), &
178               &                ldgrid = .TRUE. )
179
180            IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
181               CALL ctl_stop( 'Model not in input data' )
182               RETURN
183            ENDIF
184
185            IF ( jj == 1 ) THEN
186               ALLOCATE( clvars( inpfiles(jj)%nvar ) )
187               DO ji = 1, inpfiles(jj)%nvar
188                 clvars(ji) = inpfiles(jj)%cname(ji)
189               END DO
190            ELSE
191               DO ji = 1, inpfiles(jj)%nvar
192                  IF ( inpfiles(jj)%cname(ji) /= clvars(ji) ) THEN
193                     CALL ctl_stop( 'Feedback file variables not consistent', &
194                        &           ' with previous files for this type' )
195                  ENDIF
196               END DO
197            ENDIF
198
199            !------------------------------------------------------------------
200            !  Change longitude (-180,180)
201            !------------------------------------------------------------------
202
203            DO ji = 1, inpfiles(jj)%nobs
204
205               IF ( inpfiles(jj)%plam(ji) < -180. ) &
206                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
207
208               IF ( inpfiles(jj)%plam(ji) >  180. ) &
209                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
210
211            END DO
212
213            !------------------------------------------------------------------
214            !  Calculate the date  (change eventually)
215            !------------------------------------------------------------------
216            clrefdate=inpfiles(jj)%cdjuldref(1:8)
217            READ(clrefdate,'(I8)') irefdate(jj)
218
219            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
220            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
221               &           krefdate = irefdate(jj) )
222            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
223            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
224               &           krefdate = irefdate(jj) )
225
226            IF ( ldnightav ) THEN
227
228               IF ( lwp ) THEN
229                  WRITE(numout,*)'Resetting time of night-time averaged observations', &
230                     &             ' to the end of the day'
231               ENDIF
232
233               DO ji = 1, inpfiles(jj)%nobs
234                  !  for night-time averaged data force the time
235                  !  to be the last time-step of the day, but still within the day.
236                  IF ( inpfiles(jj)%ptim(ji) >= 0. ) THEN
237                     inpfiles(jj)%ptim(ji) = &
238                        & INT(inpfiles(jj)%ptim(ji)) + 0.9999
239                  ELSE
240                     inpfiles(jj)%ptim(ji) = &
241                        & INT(inpfiles(jj)%ptim(ji)) - 0.0001
242                  ENDIF
243               END DO
244            ENDIF
245
246            IF ( inpfiles(jj)%nobs > 0 ) THEN
247               inpfiles(jj)%iproc = -1
248               inpfiles(jj)%iobsi = -1
249               inpfiles(jj)%iobsj = -1
250            ENDIF
251            inowin = 0
252            DO ji = 1, inpfiles(jj)%nobs
253               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
254                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
255                  inowin = inowin + 1
256               ENDIF
257            END DO
258            ALLOCATE( zlam(inowin)  )
259            ALLOCATE( zphi(inowin)  )
260            ALLOCATE( iobsi(inowin) )
261            ALLOCATE( iobsj(inowin) )
262            ALLOCATE( iproc(inowin) )
263            inowin = 0
264            DO ji = 1, inpfiles(jj)%nobs
265               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
266                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
267                  inowin = inowin + 1
268                  zlam(inowin) = inpfiles(jj)%plam(ji)
269                  zphi(inowin) = inpfiles(jj)%pphi(ji)
270               ENDIF
271            END DO
272
273            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
274
275            inowin = 0
276            DO ji = 1, inpfiles(jj)%nobs
277               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
278                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
279                  inowin = inowin + 1
280                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
281                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
282                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
283               ENDIF
284            END DO
285            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
286
287            DO ji = 1, inpfiles(jj)%nobs
288               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
289                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
290                  IF ( nproc == 0 ) THEN
291                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
292                  ELSE
293                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
294                  ENDIF
295                  llvalprof = .FALSE.
296                  IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
297                     & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
298                     iobs = iobs + 1
299                  ENDIF
300               ENDIF
301            END DO
302
303         ENDIF
304
305      END DO surf_files
306
307      !-----------------------------------------------------------------------
308      ! Get the time ordered indices of the input data
309      !-----------------------------------------------------------------------
310
311      !---------------------------------------------------------------------
312      !  Loop over input data files to count total number of profiles
313      !---------------------------------------------------------------------
314      iobstot = 0
315      DO jj = 1, inobf
316         DO ji = 1, inpfiles(jj)%nobs
317            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
318               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
319               iobstot = iobstot + 1
320            ENDIF
321         END DO
322      END DO
323
324      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
325         &      isurfidx(iobstot), zdat(iobstot) )
326      jk = 0
327      DO jj = 1, inobf
328         DO ji = 1, inpfiles(jj)%nobs
329            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
330               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
331               jk = jk + 1
332               ifileidx(jk) = jj
333               isurfidx(jk) = ji
334               zdat(jk)     = inpfiles(jj)%ptim(ji)
335            ENDIF
336         END DO
337      END DO
338      CALL sort_dp_indx( iobstot, &
339         &               zdat,     &
340         &               iindx   )
341
342      CALL obs_surf_alloc( surfdata, iobs, kvars, kextr, kstp, jpi, jpj )
343
344      ! Read obs/positions, QC, all variable and assign to surfdata
345
346      iobs = 0
347
348      surfdata%cvars(:)  = clvars(:)
349
350      ityp   (:) = 0
351      itypmpp(:) = 0
352
353      ioserrcount = 0
354
355      DO jk = 1, iobstot
356
357         jj = ifileidx(iindx(jk))
358         ji = isurfidx(iindx(jk))
359         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
360            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
361
362            IF ( nproc == 0 ) THEN
363               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
364            ELSE
365               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
366            ENDIF
367
368            ! Set observation information
369
370            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
371               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
372
373               iobs = iobs + 1
374
375               CALL jul2greg( isec,                   &
376                  &           imin,                   &
377                  &           ihou,                   &
378                  &           iday,                   &
379                  &           imon,                   &
380                  &           iyea,                   &
381                  &           inpfiles(jj)%ptim(ji), &
382                  &           irefdate(jj) )
383
384
385               ! Surface time coordinates
386               surfdata%nyea(iobs) = iyea
387               surfdata%nmon(iobs) = imon
388               surfdata%nday(iobs) = iday
389               surfdata%nhou(iobs) = ihou
390               surfdata%nmin(iobs) = imin
391
392               ! Surface space coordinates
393               surfdata%rlam(iobs) = inpfiles(jj)%plam(ji)
394               surfdata%rphi(iobs) = inpfiles(jj)%pphi(ji)
395
396               ! Coordinate search parameters
397               surfdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
398               surfdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
399
400               ! Instrument type
401               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
402901            IF ( ios /= 0 ) THEN
403                  IF (ioserrcount == 0) THEN
404                     CALL ctl_warn ( 'Problem converting an instrument type ', &
405                        &            'to integer. Setting type to zero' )
406                  ENDIF
407                  ioserrcount = ioserrcount + 1
408                  itype = 0
409               ENDIF
410               surfdata%ntyp(iobs) = itype
411               IF ( itype < jpsurfmaxtype + 1 ) THEN
412                  ityp(itype+1) = ityp(itype+1) + 1
413               ELSE
414                  IF(lwp)WRITE(numout,*)'WARNING:Increase jpsurfmaxtype in ',&
415                     &                  cpname
416               ENDIF
417
418               ! Bookkeeping data to match observations
419               surfdata%nsidx(iobs) = iobs
420               surfdata%nsfil(iobs) = iindx(jk)
421
422               ! QC flags
423               surfdata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
424
425               ! Observed value
426               surfdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
427
428
429               ! Model and MDT is set to fbrmdi unless read from file
430               IF ( ldmod ) THEN
431                  surfdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
432                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) THEN
433                     surfdata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,2,1)
434                     surfdata%rext(iobs,2) = inpfiles(jj)%pext(1,ji,1)
435                  ENDIF
436                ELSE
437                  surfdata%rmod(iobs,1) = fbrmdi
438                  IF ( TRIM(surfdata%cvars(1)) == 'SLA' ) surfdata%rext(iobs,:) = fbrmdi
439               ENDIF
440            ENDIF
441         ENDIF
442
443      END DO
444
445      !-----------------------------------------------------------------------
446      ! Sum up over processors
447      !-----------------------------------------------------------------------
448
449      CALL obs_mpp_sum_integer( iobs, iobsmpp )
450      CALL obs_mpp_sum_integers( ityp, itypmpp, jpsurfmaxtype + 1 )
451
452      !-----------------------------------------------------------------------
453      ! Output number of observations.
454      !-----------------------------------------------------------------------
455      IF (lwp) THEN
456
457         WRITE(numout,*)
458         WRITE(numout,'(1X,A)')TRIM( surfdata%cvars(1) )//' data'
459         WRITE(numout,'(1X,A)')'--------------'
460         DO jj = 1,8
461            IF ( itypmpp(jj) > 0 ) THEN
462               WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
463            ENDIF
464         END DO
465         WRITE(numout,'(1X,A)') &
466            & '---------------------------------------------------------------'
467         WRITE(numout,'(1X,A,I8)') &
468            & 'Total data for variable '//TRIM( surfdata%cvars(1) )// &
469            & '           = ', iobsmpp
470         WRITE(numout,'(1X,A)') &
471            & '---------------------------------------------------------------'
472         WRITE(numout,*)
473
474      ENDIF
475
476      !-----------------------------------------------------------------------
477      ! Deallocate temporary data
478      !-----------------------------------------------------------------------
479      DEALLOCATE( ifileidx, isurfidx, zdat, clvars )
480
481      !-----------------------------------------------------------------------
482      ! Deallocate input data
483      !-----------------------------------------------------------------------
484      DO jj = 1, inobf
485         IF ( inpfiles(jj)%lalloc ) THEN
486            CALL dealloc_obfbdata( inpfiles(jj) )
487         ENDIF
488      END DO
489      DEALLOCATE( inpfiles )
490
491   END SUBROUTINE obs_rea_surf
492
493END MODULE obs_read_surf
Note: See TracBrowser for help on using the repository browser.