source: branches/UKMO/r6232_tracer_advection/NEMOGCM/NEMO/OOO_SRC/ooo_read.F90 @ 9295

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

Remove svn keywords

File size: 13.2 KB
Line 
1
2MODULE ooo_read
3   !!==================================================================
4   !!                    *** MODULE ooo_read ***
5   !! Read routines : I/O for offline obs_oper
6   !!==================================================================
7
8   USE mppini
9   USE lib_mpp
10   USE in_out_manager
11   USE par_kind, ONLY: lc
12   USE netcdf
13   USE oce,     ONLY: tsn, sshn
14   USE dom_oce, ONLY: nlci, nlcj, nimpp, njmpp, tmask
15   USE par_oce, ONLY: jpi, jpj, jpk
16   USE obs_fbm, ONLY: fbimdi, fbrmdi, fbsp, fbdp
17   USE ooo_data
18
19   IMPLICIT NONE
20   PRIVATE
21
22   PUBLIC ooo_rea_dri
23
24   !! $Id$
25CONTAINS
26   SUBROUTINE ooo_rea_dri(kfile)
27      !!------------------------------------------------------------------------
28      !!             *** ooo_rea_dri ***
29      !!
30      !! Purpose : To choose appropriate read method
31      !! Method  :
32      !!
33      !! Author  : A. Ryan Oct 2013
34      !!
35      !!------------------------------------------------------------------------
36      INTEGER, INTENT(IN) :: &
37              & kfile         !: File number
38      CHARACTER(len=lc) :: &
39              & cdfilename, & !: File name
40              & cmatchname    !: Match name
41      INTEGER :: &
42              & kindex       !: File index to read
43 
44      !! Filename, index and match-up kind
45      cdfilename = TRIM(ooo_files(kfile))
46      cmatchname = TRIM(cl4_vars(kfile))
47      kindex = nn_ooo_idx(kfile)
48
49      !! Update model fields
50      !! Class 4 variables: forecast, persistence,
51      !!                    nrt_analysis, best_estimate
52      !! Feedback variables: empty string
53      IF ( (TRIM(cmatchname) == 'forecast') .OR. &
54         & (TRIM(cmatchname) == 'persistence') .OR. &
55         & (TRIM(cmatchname) == 'nrt_analysis') .OR. &
56         & (TRIM(cmatchname) == 'best_estimate').OR. &
57         & (TRIM(cmatchname) == '') ) THEN
58         CALL ooo_read_file(TRIM(cdfilename), kindex)
59         CALL ooo_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
60      ELSE IF (TRIM(cmatchname) == 'climatology') THEN
61         WRITE(numout,*) 'Interpolating climatologies'
62      ELSE IF (TRIM(cmatchname) == 'altimeter') THEN
63         CALL ooo_read_altbias(TRIM(cdfilename))
64         CALL ooo_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
65      END IF
66
67   END SUBROUTINE ooo_rea_dri
68
69   SUBROUTINE ooo_read_altbias(filename)
70      !!------------------------------------------------------------------------
71      !!                      *** ooo_read_altbias ***
72      !!
73      !! Purpose : To read altimeter bias and set tn,sn to missing values
74      !! Method  : Use subdomain indices to create start and count matrices
75      !!           for netcdf read.
76      !!
77      !! Author  : A. Ryan Sept 2012
78      !!------------------------------------------------------------------------
79      CHARACTER(len=*), INTENT(IN) :: filename
80      INTEGER                      :: ncid, &
81                                    & varid,&
82                                    & istat,&
83                                    & ntimes,&
84                                    & tdim, &
85                                    & xdim, &
86                                    & ydim, &
87                                    & zdim
88      INTEGER                      :: ii, ij, ik
89      INTEGER, DIMENSION(3)        :: start_s, &
90                                    & count_s
91      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
92      REAL(fbdp)                     :: fill_val
93
94      IF (TRIM(filename) == 'nofile') THEN
95         tsn(:,:,:,:) = fbrmdi
96         sshn(:,:) = fbrmdi 
97      ELSE
98         ! Open Netcdf file to find dimension id
99         istat = nf90_open(trim(filename),nf90_nowrite,ncid)
100         istat = nf90_inq_dimid(ncid,'x',xdim)
101         istat = nf90_inq_dimid(ncid,'y',ydim)
102         istat = nf90_inq_dimid(ncid,'deptht',zdim)
103         istat = nf90_inq_dimid(ncid,'time',tdim)
104         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
105         ! Allocate temporary temperature array
106         ALLOCATE(temp_sshn(nlci,nlcj))
107         ! Create start and count arrays
108         start_s = (/ nimpp, njmpp, 1 /)
109         count_s = (/ nlci,  nlcj,  1 /)
110         
111         ! Altimeter bias
112         istat = nf90_inq_varid(ncid,'altbias',varid)
113         istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
114         istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
115         WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
116   
117         ! Initialise tsn, sshn to fbrmdi
118         tsn(:,:,:,:) = fbrmdi
119         sshn(:,:) = fbrmdi 
120
121         ! Fill sshn with altimeter bias
122         sshn(1:nlci,1:nlcj) = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
123
124         ! Remove halo from tmask, sshn to prevent double obs counting
125         IF (jpi > nlci) THEN
126             tmask(nlci+1:,:,:) = 0
127             sshn(nlci+1:,:) = 0
128         END IF
129         IF (jpj > nlcj) THEN
130             tmask(:,nlcj+1:,:) = 0
131             sshn(:,nlcj+1:) = 0
132         END IF
133
134         ! Deallocate arrays
135         DEALLOCATE(temp_sshn)
136
137         ! Close netcdf file
138         istat = nf90_close(ncid)
139      END IF
140   
141   END SUBROUTINE ooo_read_altbias
142
143   SUBROUTINE ooo_read_file(filename, ifcst)
144      !!------------------------------------------------------------------------
145      !!             *** ooo_read_file ***
146      !!
147      !! Purpose : To fill tn and sn with dailymean field from netcdf files
148      !! Method  : Use subdomain indices to create start and count matrices
149      !!           for netcdf read.
150      !!
151      !! Author  : A. Ryan Oct 2010
152      !!------------------------------------------------------------------------
153
154      INTEGER,          INTENT(IN) :: ifcst
155      CHARACTER(len=*), INTENT(IN) :: filename
156      INTEGER                      :: ncid, &
157                                    & varid,&
158                                    & istat,&
159                                    & ntimes,&
160                                    & tdim, &
161                                    & xdim, &
162                                    & ydim, &
163                                    & zdim
164      INTEGER                      :: ii, ij, ik
165      INTEGER, DIMENSION(4)        :: start_n, &
166                                    & count_n
167      INTEGER, DIMENSION(3)        :: start_s, &
168                                    & count_s
169      REAL(fbdp), DIMENSION(:,:,:),ALLOCATABLE :: temp_tn, &
170                                              & temp_sn
171      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
172      REAL(fbdp)                     :: fill_val
173
174      ! DEBUG
175      INTEGER :: istage
176
177      IF (TRIM(filename) == 'nofile') THEN
178         tsn(:,:,:,:) = fbrmdi
179         sshn(:,:) = fbrmdi 
180      ELSE
181         WRITE(numout,*) "Opening :", TRIM(filename)
182         ! Open Netcdf file to find dimension id
183         istat = nf90_open(path=TRIM(filename), mode=nf90_nowrite, ncid=ncid)
184         IF ( istat /= nf90_noerr ) THEN
185             WRITE(numout,*) "WARNING: Could not open ", trim(filename)
186             WRITE(numout,*) "ERROR: ", nf90_strerror(istat)
187         ENDIF
188         istat = nf90_inq_dimid(ncid,'x',xdim)
189         istat = nf90_inq_dimid(ncid,'y',ydim)
190         istat = nf90_inq_dimid(ncid,'deptht',zdim)
191         istat = nf90_inq_dimid(ncid,'time_counter',tdim)
192         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
193         IF (ifcst .LE. ntimes) THEN   
194            ! Allocate temporary temperature array
195            ALLOCATE(temp_tn(nlci,nlcj,jpk))
196            ALLOCATE(temp_sn(nlci,nlcj,jpk))
197            ALLOCATE(temp_sshn(nlci,nlcj))
198     
199            ! Set temp_tn, temp_sn to 0.
200            temp_tn(:,:,:) = fbrmdi
201            temp_sn(:,:,:) = fbrmdi
202            temp_sshn(:,:) = fbrmdi
203     
204            ! Create start and count arrays
205            start_n = (/ nimpp, njmpp, 1,   ifcst /)
206            count_n = (/ nlci,  nlcj,  jpk, 1     /)
207            start_s = (/ nimpp, njmpp, ifcst /)
208            count_s = (/ nlci,  nlcj,  1     /)
209             
210            ! Read information into temporary arrays
211            ! retrieve varid and read in temperature
212            istat = nf90_inq_varid(ncid,'votemper',varid)
213            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
214            istat = nf90_get_var(ncid, varid, temp_tn, start_n, count_n)
215            WHERE(temp_tn(:,:,:) == fill_val) temp_tn(:,:,:) = fbrmdi
216     
217            ! retrieve varid and read in salinity
218            istat = nf90_inq_varid(ncid,'vosaline',varid)
219            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
220            istat = nf90_get_var(ncid, varid, temp_sn, start_n, count_n)
221            WHERE(temp_sn(:,:,:) == fill_val) temp_sn(:,:,:) = fbrmdi
222     
223            ! retrieve varid and read in SSH
224            istat = nf90_inq_varid(ncid,'sossheig',varid)
225            IF (istat /= nf90_noerr) THEN
226               ! Altimeter bias
227               istat = nf90_inq_varid(ncid,'altbias',varid)
228            END IF
229     
230            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
231            istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
232            WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
233   
234            ! Initialise tsn, sshn to fbrmdi
235            tsn(:,:,:,:) = fbrmdi
236            sshn(:,:) = fbrmdi 
237
238            ! Mask out missing data index
239            tsn(1:nlci,1:nlcj,1:jpk,1) = temp_tn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
240            tsn(1:nlci,1:nlcj,1:jpk,2) = temp_sn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
241            sshn(1:nlci,1:nlcj)        = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
242
243            ! Remove halo from tmask, tsn, sshn to prevent double obs counting
244            IF (jpi > nlci) THEN
245                tmask(nlci+1:,:,:) = 0
246                tsn(nlci+1:,:,:,1) = 0
247                tsn(nlci+1:,:,:,2) = 0
248                sshn(nlci+1:,:) = 0
249            END IF
250            IF (jpj > nlcj) THEN
251                tmask(:,nlcj+1:,:) = 0
252                tsn(:,nlcj+1:,:,1) = 0
253                tsn(:,nlcj+1:,:,2) = 0
254                sshn(:,nlcj+1:) = 0
255            END IF
256
257            ! Deallocate arrays
258            DEALLOCATE(temp_tn, temp_sn, temp_sshn)
259         ELSE   
260            ! Mark all as missing data
261            tsn(:,:,:,:) = fbrmdi
262            sshn(:,:) = fbrmdi 
263         ENDIF
264         ! Close netcdf file
265         WRITE(numout,*) "Closing :", TRIM(filename)
266         istat = nf90_close(ncid)
267      END IF
268   END SUBROUTINE ooo_read_file
269
270   SUBROUTINE ooo_read_juld(filename, ifcst, julian)
271      USE calendar
272      !!--------------------------------------------------------------------
273      !!                 *** ooo_read_juld ***
274      !!
275      !!   Purpose : To read model julian day information from file
276      !!   Author  : A. Ryan Nov 2010
277      !!--------------------------------------------------------------------
278
279      !! Routine arguments
280      CHARACTER(len=*), INTENT(IN)  :: filename
281      INTEGER,          INTENT(IN)  :: ifcst
282      REAL,             INTENT(OUT) :: julian    !: Julian day
283
284      !! Local variables
285      INTEGER :: year,  &   !: Date information
286               & month, &
287               & day,   &
288               & hour,  &
289               & minute,&
290               & second
291      INTEGER :: istat, &   !: Netcdf variables
292               & ncid,  &
293               & dimid, &
294               & varid, &
295               & ntime     
296      REAL,DIMENSION(:),ALLOCATABLE :: r_sec     !: Remainder seconds
297      CHARACTER(len=120) :: time_str  !: time string
298
299      time_str=''
300      !! Read in time_counter and remainder seconds
301      istat = nf90_open(trim(filename),nf90_nowrite,ncid)
302      istat = nf90_inq_dimid(ncid,'time_counter',dimid)
303      IF (istat /= nf90_noerr) THEN
304          istat = nf90_inq_dimid(ncid,'time',dimid)
305      ENDIF
306      istat = nf90_inquire_dimension(ncid,dimid,len=ntime)
307      istat = nf90_inq_varid(ncid,'time_counter',varid)
308      IF (istat /= nf90_noerr) THEN
309          istat = nf90_inq_dimid(ncid,'time',dimid)
310      ENDIF
311      istat = nf90_get_att(ncid,varid,'units',time_str)
312      ALLOCATE(r_sec(ntime))
313      istat = nf90_get_var(ncid,varid, r_sec)
314      istat = nf90_close(ncid)
315
316      !! Fill yyyy-mm-dd hh:mm:ss
317      !! format(('seconds since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2))
318      100 format((14x, I4.4,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2,1x,I2.2))
319      READ( time_str, 100 ) year, month, day, hour, minute, second
320
321      CALL ymds2ju(year, month, day, r_sec(ifcst), julian)
322
323      !! To take a comment from the ymds2ju subroutine
324
325   !- In the case of the Gregorian calendar we have chosen to use
326   !- the Lilian day numbers. This is the day counter which starts
327   !- on the 15th October 1582.
328   !- This is the day at which Pope Gregory XIII introduced the
329   !- Gregorian calendar.
330   !- Compared to the true Julian calendar, which starts some
331   !- 7980 years ago, the Lilian days are smaler and are dealt with
332   !- easily on 32 bit machines. With the true Julian days you can only
333   !- the fraction of the day in the real part to a precision of
334   !- a 1/4 of a day with 32 bits.
335     
336      !! The obs operator routines prefer to calculate Julian days since
337      !! 01/01/1950 00:00:00
338      !! In order to convert to the 1950 version we must adjust by the number
339      !! of days between 15th October 1582 and 1st Jan 1950
340
341      julian = julian - 134123.
342     
343      DEALLOCATE(r_sec)
344     
345   END SUBROUTINE ooo_read_juld
346
347END MODULE ooo_read 
348
Note: See TracBrowser for help on using the repository browser.