source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/SAO_SRC/sao_read.F90 @ 4849

Last change on this file since 4849 was 4849, checked in by andrewryan, 6 years ago

moved ooo_data.F90 to sao_data.F90 along with renaming its internal subroutines appropriately

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