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/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_seaice.F90 @ 3211

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