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_sst.F90 in branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sst.F90 @ 3651

Last change on this file since 3651 was 3651, checked in by cbricaud, 11 years ago

merge dev_MERCATOR_2012_rev3555 into dev_NOC_MERCATOR_2012 ; see ticket 1020
  • Property svn:keywords set to Id
File size: 26.8 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(dp), DIMENSION(:), ALLOCATABLE :: &
113         & zdat
114      LOGICAL :: llvalprof
115      TYPE(obfbdata), POINTER, DIMENSION(:) :: &
116         & inpfiles
117      REAL(dp), DIMENSION(knumfiles) :: &
118         & djulini, &
119         & djulend
120      INTEGER :: iobs
121      INTEGER :: iobstot
122      CHARACTER(len=8) :: cl_refdate
123   
124      ! Local initialization
125      iobs = 0
126 
127      !-----------------------------------------------------------------------
128      ! Check data the model part is just with feedback data files
129      !-----------------------------------------------------------------------
130      IF ( ldmod .AND. ( kformat /= 0 ) ) THEN
131         CALL ctl_stop( 'Model can only be read from feedback data' )
132         RETURN
133      ENDIF
134
135      !-----------------------------------------------------------------------
136      ! Count the number of files needed and allocate the obfbdata type
137      !-----------------------------------------------------------------------
138     
139      inobf = knumfiles
140     
141      ALLOCATE( inpfiles(inobf) )
142
143      sst_files : DO jj = 1, inobf
144
145         CALL init_obfbdata( inpfiles(jj) )
146         
147         !---------------------------------------------------------------------
148         ! Prints
149         !---------------------------------------------------------------------
150         IF(lwp) THEN
151            WRITE(numout,*)
152            WRITE(numout,*) ' obs_rea_sst : 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 in read_proffile
180            !------------------------------------------------------------------
181           
182            iflag = nf90_close( i_file_id )
183
184            !------------------------------------------------------------------
185            !  Read the profile file into inpfiles
186            !------------------------------------------------------------------
187            IF ( kformat == 0 ) THEN
188               IF(lwp) THEN
189                  WRITE(numout,*)
190                  WRITE(numout,*)'Reading from feedback file :', &
191                     &           TRIM( cfilenames(jj) )
192               ENDIF
193               CALL read_obfbdata( TRIM( cfilenames(jj) ), inpfiles(jj), &
194                  &                ldgrid = .TRUE. )
195               IF ( ldmod .AND. ( inpfiles(jj)%nadd == 0 ) ) THEN
196                  CALL ctl_stop( 'Model not in input data' )
197                  RETURN
198               ENDIF
199            ELSEIF ( kformat == 1) THEN
200               CALL read_ghrsst( TRIM( cfilenames(jj) ), inpfiles(jj), &
201               &                 numout, lwp, .TRUE. )
202            ELSE
203               CALL ctl_stop( 'File format unknown' )
204            ENDIF
205
206            !------------------------------------------------------------------
207            !  Change longitude (-180,180)
208            !------------------------------------------------------------------
209
210            DO ji = 1, inpfiles(jj)%nobs 
211
212               IF ( inpfiles(jj)%plam(ji) < -180. ) &
213                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) + 360.
214
215               IF ( inpfiles(jj)%plam(ji) >  180. ) &
216                  &   inpfiles(jj)%plam(ji) = inpfiles(jj)%plam(ji) - 360.
217
218            END DO
219
220            !------------------------------------------------------------------
221            !  Calculate the date  (change eventually)
222            !------------------------------------------------------------------
223            cl_refdate=inpfiles(jj)%cdjuldref(1:8)
224            READ(cl_refdate,'(I8)') irefdate(jj)
225           
226            CALL ddatetoymdhms( ddobsini, iyea, imon, iday, ihou, imin, isec )
227            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulini(jj), &
228               &           krefdate = irefdate(jj) )
229            CALL ddatetoymdhms( ddobsend, iyea, imon, iday, ihou, imin, isec )
230            CALL greg2jul( isec, imin, ihou, iday, imon, iyea, djulend(jj), &
231               &           krefdate = irefdate(jj) )
232            IF ( inpfiles(jj)%nobs > 0 ) THEN
233               inpfiles(jj)%iproc = -1
234               inpfiles(jj)%iobsi = -1
235               inpfiles(jj)%iobsj = -1
236            ENDIF
237            inowin = 0
238            DO ji = 1, inpfiles(jj)%nobs
239               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
240                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
241                  inowin = inowin + 1
242               ENDIF
243            END DO
244            ALLOCATE( zlam(inowin)  )
245            ALLOCATE( zphi(inowin)  )
246            ALLOCATE( iobsi(inowin) )
247            ALLOCATE( iobsj(inowin) )
248            ALLOCATE( iproc(inowin) )
249            inowin = 0
250            DO ji = 1, inpfiles(jj)%nobs
251               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
252                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
253                  inowin = inowin + 1
254                  zlam(inowin) = inpfiles(jj)%plam(ji)
255                  zphi(inowin) = inpfiles(jj)%pphi(ji)
256               ENDIF
257            END DO
258
259            CALL obs_grid_search( inowin, zlam, zphi, iobsi, iobsj, iproc, 'T' )
260
261            inowin = 0
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                  inowin = inowin + 1
266                  inpfiles(jj)%iproc(ji,1) = iproc(inowin)
267                  inpfiles(jj)%iobsi(ji,1) = iobsi(inowin)
268                  inpfiles(jj)%iobsj(ji,1) = iobsj(inowin)
269               ENDIF
270            END DO
271            DEALLOCATE( zlam, zphi, iobsi, iobsj, iproc )
272
273            DO ji = 1, inpfiles(jj)%nobs
274               IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
275                  & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
276                  IF ( nproc == 0 ) THEN
277                     IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
278                  ELSE
279                     IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
280                  ENDIF
281                  llvalprof = .FALSE.
282                  IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
283                     & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
284                     iobs = iobs + 1
285                  ENDIF
286               ENDIF
287            END DO
288
289         ENDIF
290
291      END DO sst_files
292
293      !-----------------------------------------------------------------------
294      ! Get the time ordered indices of the input data
295      !-----------------------------------------------------------------------
296
297      !---------------------------------------------------------------------
298      !  Loop over input data files to count total number of profiles
299      !---------------------------------------------------------------------
300      iobstot = 0
301      DO jj = 1, inobf
302         DO ji = 1, inpfiles(jj)%nobs
303            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
304               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
305               iobstot = iobstot + 1
306            ENDIF
307         END DO
308      END DO
309
310      ALLOCATE( iindx(iobstot), ifileidx(iobstot), &
311         &      isstidx(iobstot), zdat(iobstot) )
312      jk = 0
313      DO jj = 1, inobf
314         DO ji = 1, inpfiles(jj)%nobs
315            IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND. &
316               & ( inpfiles(jj)%ptim(ji) <= djulend(jj) )       ) THEN
317               jk = jk + 1
318               ifileidx(jk) = jj
319               isstidx(jk) = ji
320               zdat(jk)     = inpfiles(jj)%ptim(ji)
321            ENDIF
322         END DO
323      END DO
324      CALL sort_dp_indx( iobstot, &
325         &               zdat,     &
326         &               iindx   )
327     
328      CALL obs_surf_alloc( sstdata, iobs, kvars, kextr, kstp, jpi, jpj )
329     
330      ! * Read obs/positions, QC, all variable and assign to sstdata
331 
332      iobs = 0
333
334      ityp   (:) = 0
335      itypmpp(:) = 0
336     
337      DO jk = 1, iobstot
338         
339         jj = ifileidx(iindx(jk))
340         ji = isstidx(iindx(jk))
341         IF ( ( inpfiles(jj)%ptim(ji) >  djulini(jj) ) .AND.  &
342            & ( inpfiles(jj)%ptim(ji) <= djulend(jj) ) ) THEN
343           
344            IF ( nproc == 0 ) THEN
345               IF ( inpfiles(jj)%iproc(ji,1) >  nproc ) CYCLE
346            ELSE
347               IF ( inpfiles(jj)%iproc(ji,1) /= nproc ) CYCLE
348            ENDIF
349           
350            ! Set observation information
351           
352            IF ( ( inpfiles(jj)%ivlqc(1,ji,1) == 1 ) .OR. &
353               & ( inpfiles(jj)%ivlqc(1,ji,1) == 2 ) ) THEN
354
355               iobs = iobs + 1
356
357               CALL jul2greg( isec,                   &
358                  &           imin,                   &
359                  &           ihou,                   &
360                  &           iday,                   &
361                  &           imon,                   &
362                  &           iyea,                   &
363                  &           inpfiles(jj)%ptim(ji), &
364                  &           irefdate(jj) )
365
366
367               ! SST time coordinates
368               sstdata%nyea(iobs) = iyea
369               sstdata%nmon(iobs) = imon
370               sstdata%nday(iobs) = iday
371               sstdata%nhou(iobs) = ihou
372               sstdata%nmin(iobs) = imin
373               
374               ! SST space coordinates
375               sstdata%rlam(iobs) = inpfiles(jj)%plam(ji)
376               sstdata%rphi(iobs) = inpfiles(jj)%pphi(ji)
377
378               ! Coordinate search parameters
379               sstdata%mi  (iobs) = inpfiles(jj)%iobsi(ji,1)
380               sstdata%mj  (iobs) = inpfiles(jj)%iobsj(ji,1)
381               
382               ! Instrument type
383               READ( inpfiles(jj)%cdtyp(ji), '(I4)' ) itype
384               sstdata%ntyp(iobs) = itype
385               IF ( itype < isstmaxtype + 1 ) THEN
386                  ityp(itype+1) = ityp(itype+1) + 1
387               ELSE
388                  IF(lwp)WRITE(numout,*)'WARNING:Increase isstmaxtype in ',&
389                     &                  cpname
390               ENDIF
391
392               ! Bookkeeping data to match observations
393               sstdata%nsidx(iobs) = iobs
394               sstdata%nsfil(iobs) = iindx(jk)
395
396               ! QC flags
397               sstdata%nqc(iobs) = inpfiles(jj)%ivqc(ji,1)
398
399               ! Observed value
400               sstdata%robs(iobs,1) = inpfiles(jj)%pob(1,ji,1)
401
402
403               ! Model and MDT is set to fbrmdi unless read from file
404               IF ( ldmod ) THEN
405                  sstdata%rmod(iobs,1) = inpfiles(jj)%padd(1,ji,1,1)
406               ELSE
407                  sstdata%rmod(iobs,1) = fbrmdi
408               ENDIF
409            ENDIF
410         ENDIF
411
412      END DO
413
414      !-----------------------------------------------------------------------
415      ! Sum up over processors
416      !-----------------------------------------------------------------------
417     
418      CALL obs_mpp_sum_integer( iobs, iobsmpp )
419     
420      !-----------------------------------------------------------------------
421      ! Output number of observations.
422      !-----------------------------------------------------------------------
423      IF (lwp) THEN
424
425         WRITE(numout,*)
426         WRITE(numout,'(1X,A)')'SST data types'
427         WRITE(numout,'(1X,A)')'--------------'
428         DO jj = 1,8
429            IF ( itypmpp(jj) > 0 ) THEN
430               WRITE(numout,'(1X,A4,I4,A3,I10)')'Type ', jj,' = ',itypmpp(jj)
431            ENDIF
432         END DO
433         WRITE(numout,'(1X,A50)')'--------------------------------------------------'
434         WRITE(numout,'(1X,A40,I10)')'Total                                 = ',iobsmpp
435         WRITE(numout,*)
436
437      ENDIF
438
439      !-----------------------------------------------------------------------
440      ! Deallocate temporary data
441      !-----------------------------------------------------------------------
442      DEALLOCATE( ifileidx, isstidx, zdat )
443
444      !-----------------------------------------------------------------------
445      ! Deallocate input data
446      !-----------------------------------------------------------------------
447      DO jj = 1, inobf
448         IF ( inpfiles(jj)%lalloc ) THEN
449            CALL dealloc_obfbdata( inpfiles(jj) )
450         ENDIF
451      END DO
452      DEALLOCATE( inpfiles )
453
454   END SUBROUTINE obs_rea_sst
455
456   SUBROUTINE obs_rea_sst_rey( sstname, cdsstfmt, sstdata, kvars, kextra, &
457      &                        kstp, ddobsini, ddobsend )
458      !!---------------------------------------------------------------------
459      !!
460      !!                   *** ROUTINE obs_rea_sst ***
461      !!
462      !! ** Purpose : Read from file the pseudo SST data from Reynolds
463      !!
464      !! ** Method  :
465      !!
466      !! ** Action  :
467      !!
468      !! References :
469      !!
470      !! History : 
471      !!      ! : 
472      !!----------------------------------------------------------------------
473      !! * Modules used
474      USE par_oce          ! Ocean parameters
475   
476      !! * Arguments
477      CHARACTER(len=128), INTENT(IN) :: sstname   ! Generic file name
478      CHARACTER(len=12), INTENT(IN) :: cdsstfmt   ! Format of SST files (yearly/monthly)
479      TYPE(obs_surf), INTENT(INOUT) :: sstdata    ! SST data
480      REAL(KIND=dp), INTENT(IN) :: ddobsini    ! Obs. ini time in YYYYMMDD.HHMMSS
481      REAL(KIND=dp), INTENT(IN) :: ddobsend    ! Obs. end time in YYYYMMDD.HHMMSS
482      INTEGER, INTENT(IN) :: kvars      ! Number of variables in sstdata structures
483      INTEGER, INTENT(IN) :: kextra     ! Number of extra variables in sstdata structures
484      INTEGER, INTENT(IN) :: kstp       ! Ocean time-step index
485     
486      INTEGER :: iyear
487      INTEGER :: imon
488      INTEGER :: iday
489      INTEGER :: ihour
490      INTEGER :: imin
491      INTEGER :: isec
492      INTEGER :: ihhmmss
493      INTEGER :: iyear1
494      INTEGER :: iyear2
495      INTEGER :: imon1
496      INTEGER :: imon2
497      INTEGER :: iyearf
498      INTEGER :: imonf
499      REAL(KIND=wp) :: pjulini
500      REAL(KIND=wp) :: pjulend
501      REAL(KIND=wp) :: pjulb
502      REAL(KIND=wp) :: pjule
503      REAL(KIND=wp) :: pjul
504      INTEGER :: inumsst
505      INTEGER :: itotrec
506      INTEGER :: inumobs
507      INTEGER :: irec
508      INTEGER :: ifld
509      INTEGER :: inum
510      INTEGER :: ji, jj
511      CHARACTER(len=128) :: clname
512      CHARACTER(len=4) :: cdyear
513      CHARACTER(len=2) :: cdmon
514      REAL(kind=wp), ALLOCATABLE, DIMENSION(:,:,:) :: zsstin
515
516      IF (lwp) WRITE(numout,*)'In obs_rea_sst_rey',sstname
517
518      !-----------------------------------------------------------------------
519      ! Convert observation window to julian dates.
520      !-----------------------------------------------------------------------
521      iyear1 = NINT( ddobsini / 10000 )
522      imon1 = NINT( ( ddobsini - iyear1 * 10000 ) / 100 )
523      iday = MOD( NINT( ddobsini ), 100 )
524      ihhmmss = ( ddobsini - NINT( ddobsini ) ) * 1000000
525      ihour = ihhmmss / 10000
526      imin = ( ihhmmss - ihour * 100 ) / 100
527      isec = MOD( ihhmmss, 100 )
528      CALL greg2jul ( isec, imin, ihour, iday, imon1, iyear1, pjulini )
529      IF (lwp) WRITE(numout,*)'dateini',ddobsini,iyear1,imon1,iday,ihour, &
530         & imin,isec,pjulini
531
532      iyear2 = NINT( ddobsini / 10000 )
533      imon2 = NINT( ( ddobsend - iyear2 * 10000 ) / 100 )
534      iday = MOD( NINT( ddobsend ), 100 )
535      ihhmmss = ( ddobsend - NINT( ddobsend ) ) * 1000000
536      ihour = ihhmmss / 10000
537      imin = ( ihhmmss - ihour * 100 ) / 100
538      isec = MOD( ihhmmss, 100 )
539      CALL greg2jul ( isec, imin, ihour, iday, imon2, iyear2, pjulend )
540      IF (lwp) WRITE(numout,*)'dateend',ddobsend,iyear2,imon2,iday,ihour, &
541         & imin,isec,pjulend
542
543      itotrec = NINT( pjulend - pjulini ) 
544      ALLOCATE( &
545         & zsstin( jpi, jpj, itotrec) &
546         & )
547     
548      pjul = pjulini + 1
549     
550      iyearf = -1 
551      imonf = -1
552
553      IF ( TRIM(cdsstfmt) == 'yearly' ) THEN
554
555         DO
556         
557            CALL jul2greg( isec, imin, ihour, iday, imon, iyear, &
558               &           pjul, 19500101 )
559            !
560            IF ( iyear /= iyearf ) THEN
561               
562               CALL greg2jul ( 0, 0, 0, 1, 1, iyear, pjulb )
563               
564               IF ( iyearf /= -1 ) THEN
565                 
566                  CALL iom_close ( inumsst )     
567                 
568               ENDIF
569               
570               clname = sstname
571               jj = INDEX( clname, 'YYYY' )
572               
573               IF ( jj == 0 ) THEN
574                 
575                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
576                  &           'Error processing filename ' // TRIM(sstname) )
577                 
578               ENDIF
579               
580               WRITE(cdyear,'(I4.4)') iyear
581               clname(jj:jj+3) = cdyear
582               IF(lwp) WRITE(numout,*)'Reading from Reynolds SST file : ',&
583                  & TRIM(clname)
584               
585               inumsst = 0
586               
587               CALL iom_open ( clname, inumsst )
588               
589               IF ( inumsst == 0 ) THEN
590                 
591                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
592                     &           'Error reading ' // TRIM(clname) )
593                 
594               ENDIF
595               
596               iyearf = iyear
597               
598            ENDIF
599           
600            irec = pjul - pjulb + 1
601            ifld = pjul - pjulini
602           
603            CALL iom_get ( inumsst, jpdom_data, 'sst', zsstin(:,:,ifld), irec )
604           
605            pjul = pjul + 1
606           
607            IF ( pjul > pjulend ) EXIT
608         
609         END DO
610
611      ELSEIF ( TRIM(cdsstfmt) == 'monthly' ) THEN
612
613         DO
614         
615            CALL jul2greg( isec, imin, ihour, iday, imon, iyear, &
616               &           pjul, 19500101 )
617            !
618            IF ( iyear /= iyearf .OR. imon /= imonf ) THEN
619               
620               CALL greg2jul ( 0, 0, 0, 1, imon, iyear, pjulb )
621               
622               IF ( iyearf /= -1 .AND. imonf /= -1 ) THEN
623                 
624                  CALL iom_close ( inumsst )     
625                 
626               ENDIF
627               
628               clname = sstname
629
630               jj = INDEX( clname, 'YYYY' )
631
632               IF ( jj == 0 ) THEN
633                 
634                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
635                  &           'Error processing filename ' // TRIM(sstname) )
636                 
637               ENDIF
638               
639               WRITE(cdyear,'(I4.4)') iyear
640               clname(jj:jj+3) = cdyear
641
642               jj = INDEX( clname, 'MM' )
643               
644               IF ( jj == 0 ) THEN
645                 
646                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
647                  &           'Error processing filename ' // TRIM(sstname) )
648                 
649               ENDIF
650               
651               WRITE(cdmon,'(I2.2)') imon
652               clname(jj:jj+1) = cdmon
653
654
655               IF(lwp) WRITE(numout,*)'Reading from Reynolds SST file : ',&
656                  & TRIM(clname)
657               
658               inumsst = 0
659               
660               CALL iom_open ( clname, inumsst )
661               
662               IF ( inumsst == 0 ) THEN
663                 
664                  CALL ctl_stop( 'obs_rea_sst_rey : ', &
665                     &           'Error reading ' // TRIM(clname) )
666                 
667               ENDIF
668               
669               iyearf = iyear
670               imonf = iyear
671               
672            ENDIF
673           
674            irec = pjul - pjulb + 1
675            ifld = pjul - pjulini
676           
677            CALL iom_get ( inumsst, jpdom_data, 'sst', zsstin(:,:,ifld), irec )
678           
679            pjul = pjul + 1
680           
681            IF ( pjul > pjulend ) EXIT
682         
683         END DO
684
685      ELSE
686         
687         CALL ctl_stop('Unknown REYNOLDS sst input data file format')
688
689      ENDIF
690
691      CALL iom_close ( inumsst )     
692
693      inumobs = 0
694      DO jj = nldj, nlej
695         DO ji = nldi, nlei
696            IF ( tmask(ji,jj,1) == 1.0_wp ) inumobs = inumobs + 1
697         END DO
698      END DO
699      inumobs = inumobs * itotrec
700
701      ! Allocate obs_surf data structure for time sorted data
702         
703      CALL obs_surf_alloc( sstdata, inumobs, kvars, kextra, kstp, jpi, jpj )
704
705      pjul = pjulini + 1
706
707      inumobs = 0
708
709      DO
710
711         CALL jul2greg( isec, imin, ihour, iday, imon, iyear, &
712            &           pjul, 19500101 )
713
714         ifld = pjul - pjulini
715
716         DO jj = nldj, nlej
717            DO ji = nldi, nlei
718
719               IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
720
721                  inumobs = inumobs + 1
722                 
723                  ! Integer values
724                  IF (ln_grid_global) THEN
725                     sstdata%mi(inumobs)     = MAX(mig(ji),2)
726                     sstdata%mj(inumobs)     = MAX(mjg(jj),2)
727                  ELSE
728                     sstdata%mi(inumobs)     = MAX(ji,2)
729                     sstdata%mj(inumobs)     = MAX(jj,2)
730                  ENDIF
731                  sstdata%nsidx(inumobs)  = 0
732                  sstdata%nsfil(inumobs)  = 0
733                  sstdata%nyea(inumobs)   = iyear
734                  sstdata%nmon(inumobs)   = imon
735                  sstdata%nday(inumobs)   = iday
736                  sstdata%nhou(inumobs)   = ihour
737                  sstdata%nmin(inumobs)   = imin
738                  sstdata%mstp(inumobs)   = 0
739                  sstdata%nqc(inumobs)    = 0
740                  sstdata%ntyp(inumobs)   = 0
741         
742                  ! Real values
743                  sstdata%rlam(inumobs)   = glamt(ji,jj)
744                  sstdata%rphi(inumobs)   = gphit(ji,jj)
745                  sstdata%robs(inumobs,1) = zsstin(ji,jj,ifld)
746                  sstdata%rmod(inumobs,1) = fbrmdi
747                  sstdata%rext(inumobs,:) = fbrmdi
748
749               ENDIF
750
751            END DO
752         END DO
753
754         pjul = pjul + 1
755
756         IF ( pjul > pjulend ) EXIT
757
758      END DO
759
760   END SUBROUTINE obs_rea_sst_rey
761
762END MODULE obs_read_sst
Note: See TracBrowser for help on using the repository browser.