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_seaice.F90 in branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_seaice.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

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