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 on ticket/0681/Review – Attachment – NEMO

ticket/0681/Review: obs_read_seaice.F90

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