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 @ 4094

Last change on this file since 4094 was 4084, checked in by andrewryan, 11 years ago

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