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_logchl.F90 in branches/UKMO/dev_r4650_general_vert_coord_obsoper_logchl_std/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper_logchl_std/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_logchl.F90 @ 8954

Last change on this file since 8954 was 8954, checked in by dford, 6 years ago

Import changes from reanalysis branch.
File size: 18.9 KB
Line 
1MODULE obs_read_logchl
2   !!======================================================================
3   !!                       ***  MODULE obs_read_logchl  ***
4   !! Observation diagnostics: Read the along track logchl data from
5   !!                          GHRSST or any logchl data from feedback files
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   obs_rea_logchl : Driver for reading logchl data from the 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_logchl_io            ! I/O for logchl files
24   USE iom                      ! I/O
25   USE netcdf                   ! NetCDF library
26
27   IMPLICIT NONE
28
29   !! * Routine accessibility
30   PRIVATE
31
32   PUBLIC obs_rea_logchl      ! Read the logchl 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_logchl( kformat, &
43      &                    logchldata, knumfiles, cfilenames, &
44      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
45      &                    ldignmis, ldmod )
46      !!---------------------------------------------------------------------
47      !!
48      !!                   *** ROUTINE obs_rea_logchl ***
49      !!
50      !! ** Purpose : Read from file the logchl data
51      !!
52      !! ** Method  : Depending on kformat either old or new style
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: New-style feedback
66      !                    ! 1: Old-style feedback
67      TYPE(obs_surf), INTENT(INOUT) :: &
68         & logchldata     ! logchl 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 logchldata
72      INTEGER, INTENT(IN) :: kextr    ! Number of extra fields for each var in logchldata
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_logchl'
81      INTEGER :: ji
82      INTEGER :: jj
83      INTEGER :: jk
84      INTEGER :: jadd
85      INTEGER :: iflag
86      INTEGER :: inobf
87      INTEGER :: i_file_id
88      INTEGER :: inowin
89      INTEGER :: iyea
90      INTEGER :: imon
91      INTEGER :: iday
92      INTEGER :: ihou
93      INTEGER :: imin
94      INTEGER :: isec
95      INTEGER, DIMENSION(knumfiles) :: &
96         & irefdate
97      INTEGER :: iobsmpp
98      INTEGER, PARAMETER :: ilogchlmaxtype = 1024
99      INTEGER, DIMENSION(0:ilogchlmaxtype) :: &
100         & ityp, &
101         & itypmpp
102      INTEGER, DIMENSION(:), ALLOCATABLE :: &
103         & iobsi,    &
104         & iobsj,    &
105         & iproc,    &
106         & iindx,    &
107         & ifileidx, &
108         & ilogchlidx
109      INTEGER :: itype
110      REAL(wp), DIMENSION(:), ALLOCATABLE :: &
111         & zphi, &
112         & zlam
113      real(wp), DIMENSION(:), ALLOCATABLE :: &
114         & zdat
115      LOGICAL :: llvalprof
116      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
117         & inpfiles
118      LOGICAL :: ll_get_std   ! Logical for getting STD
119      real(wp), DIMENSION(knumfiles) :: &
120         & djulini, &
121         & djulend
122      INTEGER :: iobs
123      INTEGER :: iobstot
124      INTEGER :: iaddref_std
125      INTEGER :: ios
126      INTEGER :: ioserrcount
127      CHARACTER(len=8) :: cl_refdate
128   
129      ! Local initialization
130      iobs = 0
131 
132      !-----------------------------------------------------------------------
133      ! Check data the model part is just with feedback data files
134      !-----------------------------------------------------------------------
135      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
136         CALL ctl_stop( 'Model can only be read from feedback data' )
137         RETURN
138      ENDIF
139
140      !-----------------------------------------------------------------------
141      ! Count the number of files needed and allocate the obfbdata type
142      !-----------------------------------------------------------------------
143     
144      inobf = knumfiles
145     
146      ALLOCATE( inpfiles(inobf) )
147
148      logchl_files : DO jj = 1, inobf
149         
150         !---------------------------------------------------------------------
151         ! Prints
152         !---------------------------------------------------------------------
153         IF(lwp) THEN
154            WRITE(numout,*)
155            WRITE(numout,*) ' obs_rea_logchl : Reading from file = ', &
156               & TRIM( TRIM( cfilenames(jj) ) )
157            WRITE(numout,*) ' ~~~~~~~~~~~~~~'
158            WRITE(numout,*)
159         ENDIF
160
161         !---------------------------------------------------------------------
162         !  Initialization: Open file and get dimensions only
163         !---------------------------------------------------------------------
164         
165         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
166            &                      i_file_id )
167         
168         IF ( iflag /= nf90_noerr ) THEN
169
170            IF ( ldignmis ) THEN
171               inpfiles(jj)%nobs = 0
172               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
173                  &           ' not found' )
174            ELSE
175               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
176                  &           ' not found' )
177            ENDIF
178
179         ELSE 
180           
181            !------------------------------------------------------------------
182            !  Close the file since it is opened elsewhere
183            !------------------------------------------------------------------
184           
185            iflag = nf90_close( i_file_id )
186
187            !------------------------------------------------------------------
188            !  Read the file into inpfiles
189            !------------------------------------------------------------------
190            IF ( kformat == 0 ) THEN
191               CALL init_obfbdata( inpfiles(jj) )
192               IF(lwp) THEN
193                  WRITE(numout,*)
194                  WRITE(numout,*)'Reading from feedback file :', &
195                     &           TRIM( cfilenames(jj) )
196               ENDIF
197               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
198                  &                ldgrid = .TRUE. )
199               IF ( ldmod .AND. ( ( inpfiles(jj)%nadd == 0 ) .OR.&
200                  &               ( inpfiles(jj)%next < 2 ) ) ) THEN
201                  CALL ctl_stop( 'Model not in input data' )
202                  RETURN
203               ENDIF
204            ELSEIF ( kformat == 1) THEN
205               CALL read_logchl( TRIM( cfilenames(jj) ), inpfiles(jj), &
206               &                 numout, lwp, .TRUE. )
207            ELSE
208               CALL ctl_stop( 'File format unknown' )
209            ENDIF
210
211            !------------------------------------------------------------------
212            ! Find the references to additional entries, if any
213            !------------------------------------------------------------------
214            ! STD
215            iaddref_std = 0 
216            ll_get_std = .FALSE.
217            DO jadd = 1,inpfiles(jj) % nadd 
218               IF (  TRIM(inpfiles(jj) % caddname(jadd)) == "STD" ) THEN
219                  iaddref_std = jadd
220                  ll_get_std = .TRUE. 
221                  EXIT
222               ENDIF
223            END DO
224
225            !------------------------------------------------------------------
226            !  Change longitude (-180,180)
227            !------------------------------------------------------------------
228
229            DO ji = 1, inpfiles(jj)%nobs 
230
231               IF ( inpfiles(jj)%plam(ji) < -180. ) &
232                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
233
234               IF ( inpfiles(jj)%plam(ji) >  180. ) &
235                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
236
237            END DO
238
239            !------------------------------------------------------------------
240            !  Calculate the date  (change eventually)
241            !------------------------------------------------------------------
242            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
243            READ(cl_refdate,'(I8)') irefdate(jj)
244           
245            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
246            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
247               &           krefdate = irefdate(jj) )
248            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
249            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
250               &           krefdate = irefdate(jj) )
251            IF ( inpfiles(jj)%nobs > 0 ) THEN
252               inpfiles(jj)%iproc = -1
253               inpfiles(jj)%iobsi = -1
254               inpfiles(jj)%iobsj = -1
255            ENDIF
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               ENDIF
262            END DO
263            ALLOCATE( zlam(inowin)  )
264            ALLOCATE( zphi(inowin)  )
265            ALLOCATE( iobsi(inowin) )
266            ALLOCATE( iobsj(inowin) )
267            ALLOCATE( iproc(inowin) )
268            inowin = 0
269            DO ji = 1, inpfiles(jj)%nobs
270               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
271                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
272                  inowin = inowin + 1
273                  zlam(inowin) = inpfiles(jj)%plam(ji)
274                  zphi(inowin) = inpfiles(jj)%pphi(ji)
275               ENDIF
276            END DO
277
278            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
279
280            inowin = 0
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                  inowin = inowin + 1
285                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
286                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
287                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
288               ENDIF
289            END DO
290            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
291
292            DO ji = 1, inpfiles(jj)%nobs
293               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
294                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
295                  IF ( nproc == 0 ) THEN
296                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
297                  ELSE
298                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
299                  ENDIF
300                  llvalprof = .FALSE.
301                  IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
302                     & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
303                     iobs = iobs + 1
304                  ENDIF
305               ENDIF
306            END DO
307
308         ENDIF
309
310      END DO logchl_files
311
312      !-----------------------------------------------------------------------
313      ! Get the time ordered indices of the input data
314      !-----------------------------------------------------------------------
315
316      !---------------------------------------------------------------------
317      !  Loop over input data files to count total number of obs
318      !---------------------------------------------------------------------
319      iobstot = 0
320      DO jj = 1, inobf
321         DO ji = 1, inpfiles(jj)%nobs
322            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
323               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
324               iobstot = iobstot + 1
325            ENDIF
326         END DO
327      END DO
328
329      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
330         &      ilogchlidx(iobstot), zdat(iobstot) )
331      jk = 0
332      DO jj = 1, inobf
333         DO ji = 1, inpfiles(jj)%nobs
334            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
335               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
336               jk = jk + 1
337               ifileidx(jk) = jj
338               ilogchlidx(jk) = ji
339               zdat(jk)     = inpfiles(jj)%ptim(ji)
340            ENDIF
341         END DO
342      END DO
343      CALL sort_dp_indx( iobstot, &
344         &               zdat,     &
345         &               iindx   )
346     
347      CALL obs_surf_alloc( logchldata, iobs, & 
348                           kvars, kextr, kstp, jpi, jpj )
349     
350      ! * Read obs/positions, QC, all variable and assign to logchldata
351 
352      iobs = 0
353
354      ityp   (:) = 0
355      itypmpp(:) = 0
356
357      ioserrcount=0     
358
359      DO jk = 1, iobstot
360         
361         jj = ifileidx(iindx(jk))
362         ji = ilogchlidx(iindx(jk))
363         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
364            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
365           
366            IF ( nproc == 0 ) THEN
367               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
368            ELSE
369               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
370            ENDIF
371           
372            ! Set observation information
373           
374            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
375               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
376
377               iobs = iobs + 1
378
379               CALL jul2greg( isec,                   &
380                  &           imin,                   &
381                  &           ihou,                   &
382                  &           iday,                   &
383                  &           imon,                   &
384                  &           iyea,                   &
385                  &           inpfiles(jj)%ptim(ji), &
386                  &           irefdate(jj) )
387
388
389               ! logchl time coordinates
390               logchldata%nyea(iobs) = iyea
391               logchldata%nmon(iobs) = imon
392               logchldata%nday(iobs) = iday
393               logchldata%nhou(iobs) = ihou
394               logchldata%nmin(iobs) = imin
395               
396               ! logchl space coordinates
397               logchldata%rlam(iobs) = inpfiles(jj)%plam(ji)
398               logchldata%rphi(iobs) = inpfiles(jj)%pphi(ji)
399
400               ! Coordinate search parameters
401               logchldata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
402               logchldata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
403               
404               ! Instrument type
405               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
406901            IF ( ios /= 0 ) THEN
407                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' ) 
408                  ioserrcount = ioserrcount + 1
409                  itype = 0
410               ENDIF
411               logchldata%ntyp(iobs) = itype
412               IF ( itype < ilogchlmaxtype + 1 ) THEN
413                  ityp(itype+1) = ityp(itype+1) + 1
414               ELSE
415                  IF(lwp)WRITE(numout,*)'WARNING:Increase ilogchlmaxtype in ',&
416                     &                  cpname
417               ENDIF
418
419               ! Bookkeeping data to match observations
420               logchldata%nsidx(iobs) = iobs
421               logchldata%nsfil(iobs) = iindx(jk)
422
423               ! QC flags
424               logchldata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
425
426               ! Observed value
427               logchldata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
428
429
430               ! Model and MDT is set to fbrmdi unless read from file
431               IF ( ldmod ) THEN
432                  logchldata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
433               ELSE
434                  logchldata%rmod(iobs,1) = fbrmdi
435               ENDIF
436                 
437               ! Copy in STD
438               IF ( ll_get_std ) THEN
439                  logchldata%rext(iobs,1) = inpfiles(jj)%padd(1,ji,iaddref_std,1) 
440               ELSE
441                  IF ( jk == 1 )   &
442                  &   CALL ctl_warn(   &
443                  &   "No STD values for logchl observations, setting to zero.") 
444                  logchldata%rext(iobs,1) = 0. 
445               ENDIF
446            ENDIF
447         ENDIF
448
449      END DO
450
451      !-----------------------------------------------------------------------
452      ! Sum up over processors
453      !-----------------------------------------------------------------------
454     
455      CALL obs_mpp_sum_integer( iobs, iobsmpp )
456     
457      !-----------------------------------------------------------------------
458      ! Output number of observations.
459      !-----------------------------------------------------------------------
460      IF (lwp) THEN
461
462         WRITE(numout,*)
463         WRITE(numout,'(1X,A)')'logchl data types'
464         WRITE(numout,'(1X,A)')'-----------------'
465         DO jj = 1,8
466            IF ( itypmpp(jj) > 0 ) THEN
467               WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
468            ENDIF
469         END DO
470         WRITE(numout,'(1X,A50)')'--------------------------------------------------'
471         WRITE(numout,'(1X,A40,I10)')'Total                                 = ',iobsmpp
472         WRITE(numout,*)
473
474      ENDIF
475
476      !-----------------------------------------------------------------------
477      ! Deallocate temporary data
478      !-----------------------------------------------------------------------
479      DEALLOCATE( ifileidx, ilogchlidx, zdat )
480
481      !-----------------------------------------------------------------------
482      ! Deallocate input data
483      !-----------------------------------------------------------------------
484      DO jj = 1, inobf
485         CALL dealloc_obfbdata( inpfiles(jj) )
486      END DO
487      DEALLOCATE( inpfiles )
488
489   END SUBROUTINE obs_rea_logchl
490
491END MODULE obs_read_logchl
492
Note: See TracBrowser for help on using the repository browser.