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

source: branches/UKMO/dev_r4650_general_vert_coord_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_logchl.F90 @ 7713

Last change on this file since 7713 was 7713, checked in by dford, 8 years ago

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