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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sst.F90 @ 2281

Last change on this file since 2281 was 2281, checked in by smasson, 14 years ago

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