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.
off_read.F90 in branches/2013/dev_r3987_UKMO4_OBS/NEMOGCM/NEMO/OOO_SRC – NEMO

source: branches/2013/dev_r3987_UKMO4_OBS/NEMOGCM/NEMO/OOO_SRC/off_read.F90 @ 4048

Last change on this file since 4048 was 4048, checked in by djlea, 9 years ago

Cleaning and debugging of the observation operator. Turn off the night time averaging of SST data by default, but add a namelist option to switch it on.

File size: 12.9 KB
Line 
1
2MODULE off_read
3   !!======================================================================
4   !!                      *** MODULE off_read ***
5   !! Read routines : I/O for offline obs_oper
6   !!======================================================================
7
8USE mppini
9USE lib_mpp
10USE in_out_manager
11USE par_kind, ONLY: lc
12USE netcdf
13USE oce,     ONLY: tsn, sshn
14USE dom_oce, ONLY: nlci, nlcj, nimpp, njmpp, tmask
15USE par_oce, ONLY: jpi, jpj, jpk
16USE obs_fbm, ONLY: fbimdi, fbrmdi, fbsp, fbdp
17
18USE off_data
19!! * Routine accessibility
20PRIVATE
21
22PUBLIC off_rea_dri
23
24CONTAINS
25   SUBROUTINE off_rea_dri(kfile)
26      IMPLICIT NONE
27      INTEGER, INTENT(IN) :: &
28              & kfile         !: File number
29      CHARACTER(len=lc) :: &
30              & cdfilename, & !: File name
31              & cmatchname    !: Match name
32      INTEGER :: &
33              & kindex       !: File index to read
34 
35      !! Filename, index and match-up kind
36      cdfilename = TRIM(off_files(kfile))
37      cmatchname = TRIM(cl4_vars(kfile))
38      kindex = nn_off_idx(kfile)
39
40      !! Update model fields
41      !! Class 4 variables: forecast, persistence,
42      !!                    nrt_analysis, best_estimate
43      !! Feedback variables: empty string
44      IF ( (TRIM(cmatchname) == 'forecast') .OR. &
45         & (TRIM(cmatchname) == 'persistence') .OR. &
46         & (TRIM(cmatchname) == 'nrt_analysis') .OR. &
47         & (TRIM(cmatchname) == 'best_estimate').OR. &
48         & (TRIM(cmatchname) == '') ) THEN
49         CALL off_read_file(TRIM(cdfilename), kindex)
50         CALL off_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
51      ELSE IF (TRIM(cmatchname) == 'climatology') THEN
52         WRITE(numout,*) 'Interpolating climatologies'
53      ELSE IF (TRIM(cmatchname) == 'altimeter') THEN
54         CALL off_read_altbias(TRIM(cdfilename))
55         CALL off_read_juld(TRIM(cdfilename), kindex, cl4_modjuld)
56      END IF
57
58   END SUBROUTINE off_rea_dri
59
60   SUBROUTINE off_read_altbias(filename)
61      IMPLICIT NONE
62      !------------------------------------------------------------------------
63      !             ** off_read_altbias **
64      !
65      ! Purpose : To read altimeter bias and set tn,sn to missing values
66      ! Method  : Use subdomain indices to create start and count matrices
67      !           for netcdf read.
68      !
69      ! Author  : Andrew Ryan Sept 2012
70      !
71      !------------------------------------------------------------------------
72      CHARACTER(len=*), INTENT(IN) :: filename
73      INTEGER                      :: ncid, &
74                                    & varid,&
75                                    & istat,&
76                                    & ntimes,&
77                                    & tdim, &
78                                    & xdim, &
79                                    & ydim, &
80                                    & zdim
81      INTEGER                      :: ii, ij, ik
82      INTEGER, DIMENSION(3)        :: start_s, &
83                                    & count_s
84      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
85      REAL(fbdp)                     :: fill_val
86
87      IF (TRIM(filename) == 'nofile') THEN
88         tsn(:,:,:,:) = fbrmdi
89         sshn(:,:) = fbrmdi 
90      ELSE
91         ! Open Netcdf file to find dimension id
92         istat = nf90_open(trim(filename),nf90_nowrite,ncid)
93         istat = nf90_inq_dimid(ncid,'x',xdim)
94         istat = nf90_inq_dimid(ncid,'y',ydim)
95         istat = nf90_inq_dimid(ncid,'deptht',zdim)
96         istat = nf90_inq_dimid(ncid,'time',tdim)
97         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
98         ! Allocate temporary temperature array
99         ALLOCATE(temp_sshn(nlci,nlcj))
100         ! Create start and count arrays
101         start_s = (/ nimpp, njmpp, 1 /)
102         count_s = (/ nlci,  nlcj,  1 /)
103         
104         ! Altimeter bias
105         istat = nf90_inq_varid(ncid,'altbias',varid)
106         istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
107         istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
108         WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
109   
110         ! Initialise tsn, sshn to fbrmdi
111         tsn(:,:,:,:) = fbrmdi
112         sshn(:,:) = fbrmdi 
113
114         ! Fill sshn with altimeter bias
115         sshn(1:nlci,1:nlcj) = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
116
117         ! Remove halo from tmask, sshn to prevent double obs counting
118         IF (jpi > nlci) THEN
119             tmask(nlci+1:,:,:) = 0
120             sshn(nlci+1:,:) = 0
121         END IF
122         IF (jpj > nlcj) THEN
123             tmask(:,nlcj+1:,:) = 0
124             sshn(:,nlcj+1:) = 0
125         END IF
126
127         ! Deallocate arrays
128         DEALLOCATE(temp_sshn)
129
130         ! Close netcdf file
131         istat = nf90_close(ncid)
132      END IF
133   
134   END SUBROUTINE off_read_altbias
135
136   SUBROUTINE off_read_file(filename, ifcst)
137      IMPLICIT NONE
138      !------------------------------------------------------------------------
139      !
140      !             ** off_read_file **
141      !
142      ! Purpose : To fill tn and sn with dailymean field from netcdf files
143      ! Method  : Use subdomain indices to create start and count matrices
144      !           for netcdf read.
145      !
146      ! Author  : Andrew Ryan Oct 2010
147      !
148      !------------------------------------------------------------------------
149
150      INTEGER,          INTENT(IN) :: ifcst
151      CHARACTER(len=*), INTENT(IN) :: filename
152      INTEGER                      :: ncid, &
153                                    & varid,&
154                                    & istat,&
155                                    & ntimes,&
156                                    & tdim, &
157                                    & xdim, &
158                                    & ydim, &
159                                    & zdim
160      INTEGER                      :: ii, ij, ik
161      INTEGER, DIMENSION(4)        :: start_n, &
162                                    & count_n
163      INTEGER, DIMENSION(3)        :: start_s, &
164                                    & count_s
165      REAL(fbdp), DIMENSION(:,:,:),ALLOCATABLE :: temp_tn, &
166                                              & temp_sn
167      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
168      REAL(fbdp)                     :: fill_val
169
170      ! DEBUG
171      INTEGER :: istage
172
173      IF (TRIM(filename) == 'nofile') THEN
174         tsn(:,:,:,:) = fbrmdi
175         sshn(:,:) = fbrmdi 
176      ELSE
177         WRITE(*,*) "Opening :", TRIM(filename)
178         ! Open Netcdf file to find dimension id
179         istat = nf90_open(path=TRIM(filename), mode=nf90_nowrite, ncid=ncid)
180         IF ( istat /= nf90_noerr ) THEN
181             WRITE(*,*) "WARNING: Could not open ", trim(filename)
182             WRITE(*,*) "ERROR: ", nf90_strerror(istat)
183         ENDIF
184         istat = nf90_inq_dimid(ncid,'x',xdim)
185         istat = nf90_inq_dimid(ncid,'y',ydim)
186         istat = nf90_inq_dimid(ncid,'deptht',zdim)
187         istat = nf90_inq_dimid(ncid,'time_counter',tdim)
188         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
189         IF (ifcst .LE. ntimes) THEN   
190            ! Allocate temporary temperature array
191            ALLOCATE(temp_tn(nlci,nlcj,jpk))
192            ALLOCATE(temp_sn(nlci,nlcj,jpk))
193            ALLOCATE(temp_sshn(nlci,nlcj))
194     
195            ! Set temp_tn, temp_sn to 0.
196            temp_tn(:,:,:) = fbrmdi
197            temp_sn(:,:,:) = fbrmdi
198            temp_sshn(:,:) = fbrmdi
199     
200            ! Create start and count arrays
201            start_n = (/ nimpp, njmpp, 1,   ifcst /)
202            count_n = (/ nlci,  nlcj,  jpk, 1     /)
203            start_s = (/ nimpp, njmpp, ifcst /)
204            count_s = (/ nlci,  nlcj,  1     /)
205             
206            ! Read information into temporary arrays
207            ! retrieve varid and read in temperature
208            istat = nf90_inq_varid(ncid,'votemper',varid)
209            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
210            istat = nf90_get_var(ncid, varid, temp_tn, start_n, count_n)
211            WHERE(temp_tn(:,:,:) == fill_val) temp_tn(:,:,:) = fbrmdi
212     
213            ! retrieve varid and read in salinity
214            istat = nf90_inq_varid(ncid,'vosaline',varid)
215            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
216            istat = nf90_get_var(ncid, varid, temp_sn, start_n, count_n)
217            WHERE(temp_sn(:,:,:) == fill_val) temp_sn(:,:,:) = fbrmdi
218     
219            ! retrieve varid and read in SSH
220            istat = nf90_inq_varid(ncid,'sossheig',varid)
221            IF (istat /= nf90_noerr) THEN
222               ! Altimeter bias
223               istat = nf90_inq_varid(ncid,'altbias',varid)
224            END IF
225     
226            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
227            istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
228            WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
229   
230            ! Initialise tsn, sshn to fbrmdi
231            tsn(:,:,:,:) = fbrmdi
232            sshn(:,:) = fbrmdi 
233
234            ! Mask out missing data index
235            tsn(1:nlci,1:nlcj,1:jpk,1) = temp_tn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
236            tsn(1:nlci,1:nlcj,1:jpk,2) = temp_sn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
237            sshn(1:nlci,1:nlcj)        = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
238
239            ! Remove halo from tmask, tsn, sshn to prevent double obs counting
240            IF (jpi > nlci) THEN
241                tmask(nlci+1:,:,:) = 0
242                tsn(nlci+1:,:,:,1) = 0
243                tsn(nlci+1:,:,:,2) = 0
244                sshn(nlci+1:,:) = 0
245            END IF
246            IF (jpj > nlcj) THEN
247                tmask(:,nlcj+1:,:) = 0
248                tsn(:,nlcj+1:,:,1) = 0
249                tsn(:,nlcj+1:,:,2) = 0
250                sshn(:,nlcj+1:) = 0
251            END IF
252
253            ! Deallocate arrays
254            DEALLOCATE(temp_tn, temp_sn, temp_sshn)
255         ELSE   
256            ! Mark all as missing data
257            tsn(:,:,:,:) = fbrmdi
258            sshn(:,:) = fbrmdi 
259         ENDIF
260         ! Close netcdf file
261         WRITE(*,*) "Closing :", TRIM(filename)
262         istat = nf90_close(ncid)
263      END IF
264   END SUBROUTINE off_read_file
265
266   SUBROUTINE off_read_juld(filename, ifcst, julian)
267      USE calendar
268      IMPLICIT NONE
269      !!====================================================================
270      !!
271      !!          ** off_read_juld **
272      !!
273      !!   Purpose : To read model julian day information from file
274      !!   Author  : Andrew Ryan Nov 2010
275      !!
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 off_read_juld
345
346END MODULE off_read 
347
Note: See TracBrowser for help on using the repository browser.