source: branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_seaice.F90 @ 3586

Last change on this file since 3586 was 3586, checked in by cbricaud, 8 years ago

add modification from dev_r3342_MERCATOR7_SST in dev_MERCATOR_2012_rev3555

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