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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_read_sst.F90 @ 4428

Last change on this file since 4428 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 27.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: