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

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

gross simplification of stand alone observation operator

File size: 6.5 KB
Line 
1MODULE sao_read
2   !!==================================================================
3   !!                    *** MODULE sao_read ***
4   !! Read routines : I/O for Stand Alone Observation operator
5   !!==================================================================
6   USE mppini
7   USE lib_mpp
8   USE in_out_manager
9   USE par_kind, ONLY: lc
10   USE netcdf
11   USE oce,     ONLY: tsn, sshn
12   USE dom_oce, ONLY: nlci, nlcj, nimpp, njmpp, tmask
13   USE par_oce, ONLY: jpi, jpj, jpk
14   USE obs_fbm, ONLY: fbimdi, fbrmdi, fbsp, fbdp
15   USE sao_data
16
17   IMPLICIT NONE
18   PRIVATE
19
20   PUBLIC sao_rea_dri
21
22CONTAINS
23   SUBROUTINE sao_rea_dri(kfile)
24      !!------------------------------------------------------------------------
25      !!             *** sao_rea_dri ***
26      !!
27      !! Purpose : To choose appropriate read method
28      !! Method  :
29      !!
30      !! Author  : A. Ryan Oct 2013
31      !!
32      !!------------------------------------------------------------------------
33      INTEGER, INTENT(IN) :: &
34              & kfile         !: File number
35      CHARACTER(len=lc) :: &
36              & cdfilename    !: File name
37      INTEGER :: &
38              & kindex        !: File index to read
39
40      cdfilename = TRIM(sao_files(kfile))
41      kindex = nn_sao_idx(kfile)
42      CALL sao_read_file(TRIM(cdfilename), kindex)
43
44   END SUBROUTINE sao_rea_dri
45
46   SUBROUTINE sao_read_file(filename, ifcst)
47      !!------------------------------------------------------------------------
48      !!             *** sao_read_file ***
49      !!
50      !! Purpose : To fill tn and sn with dailymean field from netcdf files
51      !! Method  : Use subdomain indices to create start and count matrices
52      !!           for netcdf read.
53      !!
54      !! Author  : A. Ryan Oct 2010
55      !!------------------------------------------------------------------------
56
57      INTEGER,          INTENT(IN) :: ifcst
58      CHARACTER(len=*), INTENT(IN) :: filename
59      INTEGER                      :: ncid, &
60                                    & varid,&
61                                    & istat,&
62                                    & ntimes,&
63                                    & tdim, &
64                                    & xdim, &
65                                    & ydim, &
66                                    & zdim
67      INTEGER                      :: ii, ij, ik
68      INTEGER, DIMENSION(4)        :: start_n, &
69                                    & count_n
70      INTEGER, DIMENSION(3)        :: start_s, &
71                                    & count_s
72      REAL(fbdp), DIMENSION(:,:,:),ALLOCATABLE :: temp_tn, &
73                                              & temp_sn
74      REAL(fbdp), DIMENSION(:,:),  ALLOCATABLE :: temp_sshn
75      REAL(fbdp)                     :: fill_val
76
77      ! DEBUG
78      INTEGER :: istage
79
80      IF (TRIM(filename) == 'nofile') THEN
81         tsn(:,:,:,:) = fbrmdi
82         sshn(:,:) = fbrmdi
83      ELSE
84         WRITE(numout,*) "Opening :", TRIM(filename)
85         ! Open Netcdf file to find dimension id
86         istat = nf90_open(path=TRIM(filename), mode=nf90_nowrite, ncid=ncid)
87         IF ( istat /= nf90_noerr ) THEN
88             WRITE(numout,*) "WARNING: Could not open ", trim(filename)
89             WRITE(numout,*) "ERROR: ", nf90_strerror(istat)
90         ENDIF
91         istat = nf90_inq_dimid(ncid,'x',xdim)
92         istat = nf90_inq_dimid(ncid,'y',ydim)
93         istat = nf90_inq_dimid(ncid,'deptht',zdim)
94         istat = nf90_inq_dimid(ncid,'time_counter',tdim)
95         istat = nf90_inquire_dimension(ncid, tdim, len=ntimes)
96         IF (ifcst .LE. ntimes) THEN
97            ! Allocate temporary temperature array
98            ALLOCATE(temp_tn(nlci,nlcj,jpk))
99            ALLOCATE(temp_sn(nlci,nlcj,jpk))
100            ALLOCATE(temp_sshn(nlci,nlcj))
101
102            ! Set temp_tn, temp_sn to 0.
103            temp_tn(:,:,:) = fbrmdi
104            temp_sn(:,:,:) = fbrmdi
105            temp_sshn(:,:) = fbrmdi
106
107            ! Create start and count arrays
108            start_n = (/ nimpp, njmpp, 1,   ifcst /)
109            count_n = (/ nlci,  nlcj,  jpk, 1     /)
110            start_s = (/ nimpp, njmpp, ifcst /)
111            count_s = (/ nlci,  nlcj,  1     /)
112
113            ! Read information into temporary arrays
114            ! retrieve varid and read in temperature
115            istat = nf90_inq_varid(ncid,'votemper',varid)
116            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
117            istat = nf90_get_var(ncid, varid, temp_tn, start_n, count_n)
118            WHERE(temp_tn(:,:,:) == fill_val) temp_tn(:,:,:) = fbrmdi
119
120            ! retrieve varid and read in salinity
121            istat = nf90_inq_varid(ncid,'vosaline',varid)
122            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
123            istat = nf90_get_var(ncid, varid, temp_sn, start_n, count_n)
124            WHERE(temp_sn(:,:,:) == fill_val) temp_sn(:,:,:) = fbrmdi
125
126            ! retrieve varid and read in SSH
127            istat = nf90_inq_varid(ncid,'sossheig',varid)
128            IF (istat /= nf90_noerr) THEN
129               ! Altimeter bias
130               istat = nf90_inq_varid(ncid,'altbias',varid)
131            END IF
132
133            istat = nf90_get_att(ncid, varid, '_FillValue', fill_val)
134            istat = nf90_get_var(ncid, varid, temp_sshn, start_s, count_s)
135            WHERE(temp_sshn(:,:) == fill_val) temp_sshn(:,:) = fbrmdi
136
137            ! Initialise tsn, sshn to fbrmdi
138            tsn(:,:,:,:) = fbrmdi
139            sshn(:,:) = fbrmdi
140
141            ! Mask out missing data index
142            tsn(1:nlci,1:nlcj,1:jpk,1) = temp_tn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
143            tsn(1:nlci,1:nlcj,1:jpk,2) = temp_sn(:,:,:) * tmask(1:nlci,1:nlcj,1:jpk)
144            sshn(1:nlci,1:nlcj)        = temp_sshn(:,:) * tmask(1:nlci,1:nlcj,1)
145
146            ! Remove halo from tmask, tsn, sshn to prevent double obs counting
147            IF (jpi > nlci) THEN
148                tmask(nlci+1:,:,:) = 0
149                tsn(nlci+1:,:,:,1) = 0
150                tsn(nlci+1:,:,:,2) = 0
151                sshn(nlci+1:,:) = 0
152            END IF
153            IF (jpj > nlcj) THEN
154                tmask(:,nlcj+1:,:) = 0
155                tsn(:,nlcj+1:,:,1) = 0
156                tsn(:,nlcj+1:,:,2) = 0
157                sshn(:,nlcj+1:) = 0
158            END IF
159
160            ! Deallocate arrays
161            DEALLOCATE(temp_tn, temp_sn, temp_sshn)
162         ELSE
163            ! Mark all as missing data
164            tsn(:,:,:,:) = fbrmdi
165            sshn(:,:) = fbrmdi
166         ENDIF
167         ! Close netcdf file
168         WRITE(numout,*) "Closing :", TRIM(filename)
169         istat = nf90_close(ncid)
170      END IF
171   END SUBROUTINE sao_read_file
172END MODULE sao_read
Note: See TracBrowser for help on using the repository browser.