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

source: branches/dev_1784_OBS/NEMO/OPA_SRC/OBS/obs_read_seaice.F90 @ 2001

Last change on this file since 2001 was 2001, checked in by djlea, 14 years ago

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