source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sst.F90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 27.2 KB
Line 
1MODULE obs_read_sst
2   !!======================================================================
3   !!                       ***  MODULE obs_read_sst  ***
4   !! Observation diagnostics: Read the SST data from the GHRSST database
5   !!                          or any SST data from feedback files
6   !!======================================================================
7
8   !!----------------------------------------------------------------------
9   !!   obs_rea_sst : Driver for reading SST data from the GHRSST/feedback
10   !!   obs_rea_sst_rey : Driver for reading SST data from Reynolds
11   !!----------------------------------------------------------------------
12
13   !! * Modules used   
14   USE par_kind                 ! Precision variables
15   USE in_out_manager           ! I/O manager
16   USE dom_oce                  ! Ocean space and time domain variables
17   USE obs_mpp                  ! MPP support routines for observation diagnostics
18   USE julian                   ! Julian date routines
19   USE obs_utils                ! Observation operator utility functions
20   USE obs_grid                 ! Grid search
21   USE obs_sort                 ! Sorting observation arrays
22   USE obs_surf_def             ! Surface observation definitions
23   USE obs_types                ! Observation type definitions
24   USE obs_sst_io               ! I/O for sst files
25   USE iom                      ! I/O of fields for Reynolds data
26   USE netcdf                   ! NetCDF library
27
28   IMPLICIT NONE
29
30   !! * Routine accessibility
31   PRIVATE
32
33   PUBLIC obs_rea_sst      ! Read the SST observations from the point data
34   PUBLIC obs_rea_sst_rey  ! Read the gridded Reynolds SST
35   
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41
42CONTAINS
43
44   SUBROUTINE obs_rea_sst( kformat, &
45      &                    sstdata, knumfiles, cfilenames, &
46      &                    kvars, kextr, kstp, ddobsini, ddobsend, &
47      &                    ldignmis, ldmod )
48      !!---------------------------------------------------------------------
49      !!
50      !!                   *** ROUTINE obs_rea_sst ***
51      !!
52      !! ** Purpose : Read from file the SST data
53      !!
54      !! ** Method  : Depending on kformat either AVISO or
55      !!              feedback data files are read
56      !!
57      !! ** Action  :
58      !!
59      !!
60      !! History : 
61      !!      ! :  2009-01 (K. Mogensen) Initial version based on old versions
62      !!----------------------------------------------------------------------
63      !! * Modules used
64
65      !! * Arguments
66      INTEGER :: kformat   ! Format of input data
67      !                    ! 0: Feedback
68      !                    ! 1: GHRSST
69      TYPE(obs_surf), INTENT(INOUT) :: sstdata   ! SST data to be read
70      INTEGER, INTENT(IN) :: knumfiles    ! Number of corio format files to read in
71      CHARACTER(LEN=128), INTENT(IN) :: cfilenames(knumfiles) ! File names to read in
72      INTEGER, INTENT(IN) :: kvars      ! Number of variables in sstdata
73      INTEGER, INTENT(IN) :: kextr      ! Number of extra fields for each var in sstdata
74      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
75      LOGICAL, INTENT(IN) :: ldignmis   ! Ignore missing files
76      LOGICAL, INTENT(IN) :: ldmod      ! Initialize model from input data
77      REAL(KIND=dp), INTENT(IN) :: ddobsini   ! Obs. ini time in YYYYMMDD.HHMMSS
78      REAL(KIND=dp), INTENT(IN) :: ddobsend   ! Obs. end time in YYYYMMDD.HHMMSS
79         
80      !! * Local declarations
81      CHARACTER(LEN=11), PARAMETER :: cpname='obs_rea_sst'
82      INTEGER :: ji
83      INTEGER :: jj
84      INTEGER :: jk
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) :: irefdate
96      INTEGER :: iobsmpp
97      INTEGER, PARAMETER :: isstmaxtype = 1024
98      INTEGER, DIMENSION(0:isstmaxtype) :: &
99         & ityp, &
100         & itypmpp
101      INTEGER, DIMENSION(:), ALLOCATABLE :: &
102         & iobsi,    &
103         & iobsj,    &
104         & iproc,    &
105         & iindx,    &
106         & ifileidx, &
107         & isstidx
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      sst_files : DO jj = 1, inobf
146
147         CALL init_obfbdata( inpfiles(jj) )
148         
149         !---------------------------------------------------------------------
150         ! Prints
151         !---------------------------------------------------------------------
152         IF(lwp) THEN
153            WRITE(numout,*)
154            WRITE(numout,*) ' obs_rea_sst : Reading from file = ', &
155               & TRIM( TRIM( cfilenames(jj) ) )
156            WRITE(numout,*) ' ~~~~~~~~~~~'
157            WRITE(numout,*)
158         ENDIF
159
160         !---------------------------------------------------------------------
161         !  Initialization: Open file and get dimensions only
162         !---------------------------------------------------------------------
163         
164         iflag = nf90_open( TRIM( TRIM( cfilenames(jj) ) ), nf90_nowrite, &
165            &                      i_file_id )
166         
167         IF ( iflag /= nf90_noerr ) THEN
168
169            IF ( ldignmis ) THEN
170               inpfiles(jj)%nobs = 0
171               CALL ctl_warn( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
172                  &           ' not found' )
173            ELSE
174               CALL ctl_stop( 'File ' // TRIM( TRIM( cfilenames(jj) ) ) // &
175                  &           ' not found' )
176            ENDIF
177
178         ELSE 
179           
180            !------------------------------------------------------------------
181            !  Close the file since it is opened in read_proffile
182            !------------------------------------------------------------------
183           
184            iflag = nf90_close( i_file_id )
185
186            !------------------------------------------------------------------
187            !  Read the profile file into inpfiles
188            !------------------------------------------------------------------
189            IF ( kformat == 0 ) THEN
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 ) ) THEN
198                  CALL ctl_stop( 'Model not in input data' )
199                  RETURN
200               ENDIF
201            ELSEIF ( kformat == 1) THEN
202               CALL read_ghrsst( 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 sst_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 profiles
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         &      isstidx(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               isstidx(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( sstdata, iobs, kvars, kextr, kstp, jpi, jpj )
331     
332      ! * Read obs/positions, QC, all variable and assign to sstdata
333 
334      iobs = 0
335
336      ityp   (:) = 0
337      itypmpp(:) = 0
338     
339      ioserrcount = 0
340     
341      DO jk = 1, iobstot
342         
343         jj = ifileidx(iindx(jk))
344         ji = isstidx(iindx(jk))
345         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
346            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
347           
348            IF ( nproc == 0 ) THEN
349               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
350            ELSE
351               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
352            ENDIF
353           
354            ! Set observation information
355           
356            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
357               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
358
359               iobs = iobs + 1
360
361               CALL jul2greg( isec,                   &
362                  &           imin,                   &
363                  &           ihou,                   &
364                  &           iday,                   &
365                  &           imon,                   &
366                  &           iyea,                   &
367                  &           inpfiles(jj)%ptim(ji), &
368                  &           irefdate(jj) )
369
370
371               ! SST time coordinates
372               sstdata%nyea(iobs) = iyea
373               sstdata%nmon(iobs) = imon
374               sstdata%nday(iobs) = iday
375               sstdata%nhou(iobs) = ihou
376               sstdata%nmin(iobs) = imin
377               
378               ! SST space coordinates
379               sstdata%rlam(iobs) = inpfiles(jj)%plam(ji)
380               sstdata%rphi(iobs) = inpfiles(jj)%pphi(ji)
381
382               ! Coordinate search parameters
383               sstdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
384               sstdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
385               
386               ! Instrument type
387               READ( inpfiles(jj)%cdtyp(ji), '(I4)', IOSTAT = ios, ERR = 901 ) itype
388901            IF ( ios /= 0 ) THEN
389                  IF (ioserrcount == 0) CALL ctl_warn ( 'Problem converting an instrument type to integer. Setting type to zero' ) 
390                  ioserrcount = ioserrcount + 1
391                  itype = 0
392               ENDIF
393               sstdata%ntyp(iobs) = itype
394               IF ( itype < isstmaxtype + 1 ) THEN
395                  ityp(itype+1) = ityp(itype+1) + 1
396               ELSE
397                  IF(lwp)WRITE(numout,*)'WARNING:Increase isstmaxtype in ',&
398                     &                  cpname
399               ENDIF
400
401               ! Bookkeeping data to match observations
402               sstdata%nsidx(iobs) = iobs
403               sstdata%nsfil(iobs) = iindx(jk)
404
405               ! QC flags
406               sstdata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
407
408               ! Observed value
409               sstdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
410
411
412               ! Model and MDT is set to fbrmdi unless read from file
413               IF ( ldmod ) THEN
414                  sstdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
415               ELSE
416                  sstdata%rmod(iobs,1) = fbrmdi
417               ENDIF
418            ENDIF
419         ENDIF
420
421      END DO
422
423      !-----------------------------------------------------------------------
424      ! Sum up over processors
425      !-----------------------------------------------------------------------
426     
427      CALL obs_mpp_sum_integer( iobs, iobsmpp )
428     
429      !-----------------------------------------------------------------------
430      ! Output number of observations.
431      !-----------------------------------------------------------------------
432      IF (lwp) THEN
433
434         WRITE(numout,*)
435         WRITE(numout,'(1X,A)')'SST data types'
436         WRITE(numout,'(1X,A)')'--------------'
437         DO jj = 1,8
438            IF ( itypmpp(jj) > 0 ) THEN
439               WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
440            ENDIF
441         END DO
442         WRITE(numout,'(1X,A50)')'--------------------------------------------------'
443         WRITE(numout,'(1X,A40,I10)')'Total                                 = ',iobsmpp
444         WRITE(numout,*)
445
446      ENDIF
447
448      !-----------------------------------------------------------------------
449      ! Deallocate temporary data
450      !-----------------------------------------------------------------------
451      DEALLOCATE( ifileidx, isstidx, zdat )
452
453      !-----------------------------------------------------------------------
454      ! Deallocate input data
455      !-----------------------------------------------------------------------
456      DO jj = 1, inobf
457         IF ( inpfiles(jj)%lalloc ) THEN
458            CALL dealloc_obfbdata( inpfiles(jj) )
459         ENDIF
460      END DO
461      DEALLOCATE( inpfiles )
462
463   END SUBROUTINE obs_rea_sst
464
465   SUBROUTINE obs_rea_sst_rey( sstname, cdsstfmt, sstdata, kvars, kextra, &
466      &                        kstp, ddobsini, ddobsend )
467      !!---------------------------------------------------------------------
468      !!
469      !!                   *** ROUTINE obs_rea_sst ***
470      !!
471      !! ** Purpose : Read from file the pseudo SST data from Reynolds
472      !!
473      !! ** Method  :
474      !!
475      !! ** Action  :
476      !!
477      !! References :
478      !!
479      !! History : 
480      !!      ! : 
481      !!----------------------------------------------------------------------
482      !! * Modules used
483      USE par_oce          ! Ocean parameters
484   
485      !! * Arguments
486      CHARACTER(len=128), INTENT(IN) :: sstname   ! Generic file name
487      CHARACTER(len=12), INTENT(IN) :: cdsstfmt   ! Format of SST files (yearly/monthly)
488      TYPE(obs_surf), INTENT(INOUT) :: sstdata    ! SST data
489      REAL(KIND=dp), INTENT(IN) :: ddobsini    ! Obs. ini time in YYYYMMDD.HHMMSS
490      REAL(KIND=dp), INTENT(IN) :: ddobsend    ! Obs. end time in YYYYMMDD.HHMMSS
491      INTEGER, INTENT(IN) :: kvars      ! Number of variables in sstdata structures
492      INTEGER, INTENT(IN) :: kextra     ! Number of extra variables in sstdata structures
493      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
494     
495      INTEGER :: iyear
496      INTEGER :: imon
497      INTEGER :: iday
498      INTEGER :: ihour
499      INTEGER :: imin
500      INTEGER :: isec
501      INTEGER :: ihhmmss
502      INTEGER :: iyear1
503      INTEGER :: iyear2
504      INTEGER :: imon1
505      INTEGER :: imon2
506      INTEGER :: iyearf
507      INTEGER :: imonf
508      REAL(KIND=wp) :: pjulini
509      REAL(KIND=wp) :: pjulend
510      REAL(KIND=wp) :: pjulb
511      REAL(KIND=wp) :: pjule
512      REAL(KIND=wp) :: pjul
513      INTEGER :: inumsst
514      INTEGER :: itotrec
515      INTEGER :: inumobs
516      INTEGER :: irec
517      INTEGER :: ifld
518      INTEGER :: inum
519      INTEGER :: ji, jj
520      CHARACTER(len=128) :: clname
521      CHARACTER(len=4) :: cdyear
522      CHARACTER(len=2) :: cdmon
523      REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: zsstin
524
525      IF (lwp) WRITE(numout,*)'In obs_rea_sst_rey',sstname
526
527      !-----------------------------------------------------------------------
528      ! Convert observation window to julian dates.
529      !-----------------------------------------------------------------------
530      iyear1 = NINT( ddobsini / 10000 )
531      imon1 = NINT( ( ddobsini - iyear1 * 10000 ) / 100 )
532      iday = MOD( NINT( ddobsini ), 100 )
533      ihhmmss = ( ddobsini - NINT( ddobsini ) ) * 1000000
534      ihour = ihhmmss / 10000
535      imin = ( ihhmmss - ihour * 100 ) / 100
536      isec = MOD( ihhmmss, 100 )
537      CALL greg2jul ( isec, imin, ihour, iday, imon1, iyear1, pjulini )
538      IF (lwp) WRITE(numout,*)'dateini',ddobsini,iyear1,imon1,iday,ihour, &
539         & imin,isec,pjulini
540
541      iyear2 = NINT( ddobsini / 10000 )
542      imon2 = NINT( ( ddobsend - iyear2 * 10000 ) / 100 )
543      iday = MOD( NINT( ddobsend ), 100 )
544      ihhmmss = ( ddobsend - NINT( ddobsend ) ) * 1000000
545      ihour = ihhmmss / 10000
546      imin = ( ihhmmss - ihour * 100 ) / 100
547      isec = MOD( ihhmmss, 100 )
548      CALL greg2jul ( isec, imin, ihour, iday, imon2, iyear2, pjulend )
549      IF (lwp) WRITE(numout,*)'dateend',ddobsend,iyear2,imon2,iday,ihour, &
550         & imin,isec,pjulend
551
552      itotrec = NINT( pjulend - pjulini ) 
553      ALLOCATE( &
554         & zsstin( jpi, jpj, itotrec) &
555         & )
556     
557      pjul = pjulini + 1
558     
559      iyearf = -1 
560      imonf = -1
561
562      IF ( TRIM(cdsstfmt) == 'yearly' ) THEN
563
564         DO
565         
566            CALL jul2greg( isec, imin, ihour, iday, imon, iyear, &
567               &           pjul, 19500101 )
568            !
569            IF ( iyear /= iyearf ) THEN
570               
571               CALL greg2jul ( 0, 0, 0, 1, 1, iyear, pjulb )
572               
573               IF ( iyearf /= -1 ) THEN
574                 
575                  CALL iom_close ( inumsst )     
576                 
577               ENDIF
578               
579               clname = sstname
580               jj = INDEX( clname, 'YYYY' )
581               
582               IF ( jj == 0 ) THEN
583                 
584                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
585                  &           'Error processing filename ' // TRIM(sstname) )
586                 
587               ENDIF
588               
589               WRITE(cdyear,'(I4.4)') iyear
590               clname(jj:jj+3) = cdyear
591               IF(lwp) WRITE(numout,*)'Reading from Reynolds SST file : ',&
592                  & TRIM(clname)
593               
594               inumsst = 0
595               
596               CALL iom_open ( clname, inumsst )
597               
598               IF ( inumsst == 0 ) THEN
599                 
600                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
601                     &           'Error reading ' // TRIM(clname) )
602                 
603               ENDIF
604               
605               iyearf = iyear
606               
607            ENDIF
608           
609            irec = pjul - pjulb + 1
610            ifld = pjul - pjulini
611           
612            CALL iom_get ( inumsst, jpdom_data, 'sst', zsstin(:,:,ifld), irec )
613           
614            pjul = pjul + 1
615           
616            IF ( pjul > pjulend ) EXIT
617         
618         END DO
619
620      ELSEIF ( TRIM(cdsstfmt) == 'monthly' ) THEN
621
622         DO
623         
624            CALL jul2greg( isec, imin, ihour, iday, imon, iyear, &
625               &           pjul, 19500101 )
626            !
627            IF ( iyear /= iyearf .OR. imon /= imonf ) THEN
628               
629               CALL greg2jul ( 0, 0, 0, 1, imon, iyear, pjulb )
630               
631               IF ( iyearf /= -1 .AND. imonf /= -1 ) THEN
632                 
633                  CALL iom_close ( inumsst )     
634                 
635               ENDIF
636               
637               clname = sstname
638
639               jj = INDEX( clname, 'YYYY' )
640
641               IF ( jj == 0 ) THEN
642                 
643                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
644                  &           'Error processing filename ' // TRIM(sstname) )
645                 
646               ENDIF
647               
648               WRITE(cdyear,'(I4.4)') iyear
649               clname(jj:jj+3) = cdyear
650
651               jj = INDEX( clname, 'MM' )
652               
653               IF ( jj == 0 ) THEN
654                 
655                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
656                  &           'Error processing filename ' // TRIM(sstname) )
657                 
658               ENDIF
659               
660               WRITE(cdmon,'(I2.2)') imon
661               clname(jj:jj+1) = cdmon
662
663
664               IF(lwp) WRITE(numout,*)'Reading from Reynolds SST file : ',&
665                  & TRIM(clname)
666               
667               inumsst = 0
668               
669               CALL iom_open ( clname, inumsst )
670               
671               IF ( inumsst == 0 ) THEN
672                 
673                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
674                     &           'Error reading ' // TRIM(clname) )
675                 
676               ENDIF
677               
678               iyearf = iyear
679               imonf = iyear
680               
681            ENDIF
682           
683            irec = pjul - pjulb + 1
684            ifld = pjul - pjulini
685           
686            CALL iom_get ( inumsst, jpdom_data, 'sst', zsstin(:,:,ifld), irec )
687           
688            pjul = pjul + 1
689           
690            IF ( pjul > pjulend ) EXIT
691         
692         END DO
693
694      ELSE
695         
696         CALL ctl_stop('Unknown REYNOLDS sst input data file format')
697
698      ENDIF
699
700      CALL iom_close ( inumsst )     
701
702      inumobs = 0
703      DO jj = nldj, nlej
704         DO ji = nldi, nlei
705            IF ( tmask(ji,jj,1) == 1.0_wp ) inumobs = inumobs + 1
706         END DO
707      END DO
708      inumobs = inumobs * itotrec
709
710      ! Allocate obs_surf data structure for time sorted data
711         
712      CALL obs_surf_alloc( sstdata, inumobs, kvars, kextra, kstp, jpi, jpj )
713
714      pjul = pjulini + 1
715
716      inumobs = 0
717
718      DO
719
720         CALL jul2greg( isec, imin, ihour, iday, imon, iyear, &
721            &           pjul, 19500101 )
722
723         ifld = pjul - pjulini
724
725         DO jj = nldj, nlej
726            DO ji = nldi, nlei
727
728               IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
729
730                  inumobs = inumobs + 1
731                 
732                  ! Integer values
733                  IF (ln_grid_global) THEN
734                     sstdata%mi(inumobs)     = MAX(mig(ji),2)
735                     sstdata%mj(inumobs)     = MAX(mjg(jj),2)
736                  ELSE
737                     sstdata%mi(inumobs)     = MAX(ji,2)
738                     sstdata%mj(inumobs)     = MAX(jj,2)
739                  ENDIF
740                  sstdata%nsidx(inumobs)  = 0
741                  sstdata%nsfil(inumobs)  = 0
742                  sstdata%nyea(inumobs)   = iyear
743                  sstdata%nmon(inumobs)   = imon
744                  sstdata%nday(inumobs)   = iday
745                  sstdata%nhou(inumobs)   = ihour
746                  sstdata%nmin(inumobs)   = imin
747                  sstdata%mstp(inumobs)   = 0
748                  sstdata%nqc(inumobs)    = 0
749                  sstdata%ntyp(inumobs)   = 0
750         
751                  ! Real values
752                  sstdata%rlam(inumobs)   = glamt(ji,jj)
753                  sstdata%rphi(inumobs)   = gphit(ji,jj)
754                  sstdata%robs(inumobs,1) = zsstin(ji,jj,ifld)
755                  sstdata%rmod(inumobs,1) = fbrmdi
756                  sstdata%rext(inumobs,:) = fbrmdi
757
758               ENDIF
759
760            END DO
761         END DO
762
763         pjul = pjul + 1
764
765         IF ( pjul > pjulend ) EXIT
766
767      END DO
768
769   END SUBROUTINE obs_rea_sst_rey
770
771END MODULE obs_read_sst
Note: See TracBrowser for help on using the repository browser.