source: branches/UKMO/r6232_obs_oper_update/NEMOGCM/NEMO/OPA_SRC/OBS/obs_sstbias.F90 @ 11202

Last change on this file since 11202 was 11202, checked in by jcastill, 15 months ago

Copy of branch branches/UKMO/dev_r5518_obs_oper_update@11130 without namelist_ref changes to allow merging with coupling and biogeochemistry branches

File size: 9.3 KB
Line 
1MODULE obs_sstbias
2   !!======================================================================
3   !!                       ***  MODULE obs_sstbias  ***
4   !! Observation diagnostics: Read the bias for SST data
5   !!======================================================================
6   !!----------------------------------------------------------------------
7   !!   obs_app_sstbias : Driver for reading and applying the SST bias
8   !!----------------------------------------------------------------------
9   !! * Modules used   
10   USE par_kind, ONLY : &       ! Precision variables
11      & wp, &
12      & dp, &
13      & sp
14   USE par_oce, ONLY : &        ! Domain parameters
15      & jpi, &
16      & jpj, &
17      & jpim1
18   USE in_out_manager, ONLY : & ! I/O manager
19      & lwp,    &
20      & numout 
21   USE obs_surf_def             ! Surface observation definitions
22   USE dom_oce, ONLY : &        ! Domain variables
23      & tmask, &
24      & gphit, &
25      & glamt
26   USE obs_inter_h2d
27   USE obs_utils               ! Various observation tools
28   USE obs_inter_sup
29   IMPLICIT NONE
30   !! * Routine accessibility
31   PRIVATE
32   PUBLIC obs_app_sstbias     ! Read the altimeter bias
33CONTAINS
34   SUBROUTINE obs_app_sstbias( sstdata, k2dint, knumtypes, &
35                               cl_bias_files )
36      !!---------------------------------------------------------------------
37      !!
38      !!                   *** ROUTINE obs_app_sstbias ***
39      !!
40      !! ** Purpose : Read SST bias data from files and apply correction to
41      !!             observations
42      !!
43      !! ** Method  :
44      !!
45      !! ** Action  :
46      !!
47      !! References :
48      !!
49      !! History :
50      !!      ! :  2014-08 (J. While) Bias correction code for SST obs,
51      !!      !                       based on obs_rea_altbias
52      !!----------------------------------------------------------------------
53      !! * Modules used
54      USE iom
55      USE netcdf
56
57      !! * Arguments
58      TYPE(obs_surf), INTENT(INOUT) :: &
59         & sstdata            ! SST data
60      INTEGER, INTENT(IN) :: k2dint
61      INTEGER, INTENT(IN) :: &
62         & knumtypes          ! Number of bias types to read in
63      CHARACTER(LEN=128), DIMENSION(knumtypes), INTENT(IN) :: &
64         & cl_bias_files      ! List of files to read
65
66      !! * Local declarations
67      INTEGER :: jslano       ! Data set loop variable
68      INTEGER :: jobs         ! Obs loop variable
69      INTEGER :: jpisstbias   ! Number of grid point in latitude for the bias
70      INTEGER :: jpjsstbias   ! Number of grid point in longitude for the bias
71      INTEGER :: jt
72      INTEGER :: i_nx_id      ! Index to read the NetCDF file
73      INTEGER :: i_ny_id      !
74      INTEGER :: i_file_id    !
75      INTEGER :: i_var_id
76      INTEGER, DIMENSION(knumtypes) :: &
77         & ibiastypes         ! Array of the bias types in each file
78      REAL(wp), DIMENSION(jpi,jpj,knumtypes) :: & 
79         & z_sstbias          ! Array to store the SST bias values
80      REAL(wp), DIMENSION(jpi,jpj) :: & 
81         & z_sstbias_2d       ! Array to store the SST bias values   
82      REAL(wp), DIMENSION(1) :: &
83         & zext, &
84         & zobsmask
85      REAL(wp), DIMENSION(2,2,1) :: &
86         & zweig
87      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
88         & zmask, &
89         & zglam, &
90         & zgphi
91      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
92         & zmask_tmp, &
93         & zglam_tmp, &
94         & zgphi_tmp   
95      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zbias   
96      REAL(wp) :: zlam
97      REAL(wp) :: zphi
98      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
99         & igrdi, &
100         & igrdj
101      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
102         & igrdi_tmp, &
103         & igrdj_tmp   
104      INTEGER :: numsstbias
105      INTEGER(KIND=NF90_INT) :: ifile_source
106     
107      INTEGER :: incfile
108      INTEGER :: jtype
109      INTEGER :: iret 
110      INTEGER :: inumtype
111
112      IF ( lwp ) THEN
113         WRITE(numout,*) 
114         WRITE(numout,*) 'obs_app_sstbias : '
115         WRITE(numout,*) '----------------- '
116         WRITE(numout,*) 'Read SST bias '
117      ENDIF
118
119      ! Open and read the SST bias files for each bias type
120      z_sstbias(:,:,:) = 0.0_wp
121
122      DO jtype = 1, knumtypes
123     
124         numsstbias=0
125
126         IF ( lwp ) WRITE(numout,*) 'Opening ',cl_bias_files(jtype)
127         CALL iom_open( cl_bias_files(jtype), numsstbias, ldstop=.FALSE. )
128
129         IF (numsstbias > 0) THEN
130     
131            !Read the bias type from the file
132            !No IOM get attribute command at time of writing,
133            !so have to use NETCDF
134            !routines directly - should be upgraded in the future
135            iret=NF90_OPEN(TRIM(cl_bias_files(jtype)), NF90_NOWRITE, incfile)
136            iret=NF90_GET_ATT( incfile, NF90_GLOBAL, "SST_source", &
137                               ifile_source )
138            ibiastypes(jtype) = ifile_source                 
139            iret=NF90_CLOSE(incfile)       
140           
141            IF ( iret /= 0  ) THEN
142               CALL ctl_stop( 'obs_app_sstbias : Cannot read bias type from file '// &
143                  &           TRIM( cl_bias_files(jtype) ) )
144            ENDIF
145
146            ! Get the SST bias data
147            CALL iom_get( numsstbias, jpdom_data, 'tn', z_sstbias_2d(:,:), 1 )
148            z_sstbias(:,:,jtype) = z_sstbias_2d(:,:)       
149            ! Close the file
150            CALL iom_close(numsstbias)
151            IF ( lwp ) WRITE(numout,*) 'Read in bias for type: ',ibiastypes(jtype)
152         ELSE
153            CALL ctl_stop('obs_read_sstbias: File '// & 
154               &          TRIM( cl_bias_files(jtype) )//' Not found')
155         ENDIF
156
157      END DO
158           
159      ! Interpolate the bias from the model grid to the observation points
160      ALLOCATE( &
161         & igrdi(2,2,sstdata%nsurf), &
162         & igrdj(2,2,sstdata%nsurf), &
163         & zglam(2,2,sstdata%nsurf), &
164         & zgphi(2,2,sstdata%nsurf), &
165         & zmask(2,2,sstdata%nsurf)  )
166       
167      DO jobs = 1, sstdata%nsurf 
168         igrdi(1,1,jobs) = sstdata%mi(jobs)-1
169         igrdj(1,1,jobs) = sstdata%mj(jobs)-1
170         igrdi(1,2,jobs) = sstdata%mi(jobs)-1
171         igrdj(1,2,jobs) = sstdata%mj(jobs)
172         igrdi(2,1,jobs) = sstdata%mi(jobs)
173         igrdj(2,1,jobs) = sstdata%mj(jobs)-1
174         igrdi(2,2,jobs) = sstdata%mi(jobs)
175         igrdj(2,2,jobs) = sstdata%mj(jobs)
176      END DO
177
178      CALL obs_int_comm_2d( 2, 2, sstdata%nsurf, jpi, jpj, &
179         &                  igrdi, igrdj, glamt, zglam )
180      CALL obs_int_comm_2d( 2, 2, sstdata%nsurf, jpi, jpj, &
181         &                  igrdi, igrdj, gphit, zgphi )
182      CALL obs_int_comm_2d( 2, 2, sstdata%nsurf, jpi, jpj, &
183         &                  igrdi, igrdj, tmask(:,:,1), zmask )
184
185      DO jtype = 1, knumtypes
186         
187         !Find the number observations of this type
188         !and allocate tempory arrays
189         inumtype = COUNT( sstdata%ntyp(:) == ibiastypes(jtype) )
190
191         ALLOCATE( &
192            & igrdi_tmp(2,2,inumtype), &
193            & igrdj_tmp(2,2,inumtype), &
194            & zglam_tmp(2,2,inumtype), &
195            & zgphi_tmp(2,2,inumtype), &
196            & zmask_tmp(2,2,inumtype), &
197            & zbias( 2,2,inumtype ) )
198
199         IF ( lwp ) WRITE(numout,*) 'Number of obs of type: ',ibiastypes(jtype),' is: ',inumtype
200         
201         jt=1
202         DO jobs = 1, sstdata%nsurf 
203            IF ( sstdata%ntyp(jobs) == ibiastypes(jtype) ) THEN
204
205               igrdi_tmp(:,:,jt) = igrdi(:,:,jobs) 
206               igrdj_tmp(:,:,jt) = igrdj(:,:,jobs)
207               zglam_tmp(:,:,jt) = zglam(:,:,jobs)
208               zgphi_tmp(:,:,jt) = zgphi(:,:,jobs)
209               zmask_tmp(:,:,jt) = zmask(:,:,jobs)
210
211               jt = jt +1
212
213            ENDIF
214         END DO
215                         
216         CALL obs_int_comm_2d( 2, 2, inumtype, jpi, jpj,       &
217               &           igrdi_tmp(:,:,:), igrdj_tmp(:,:,:), &
218               &           z_sstbias(:,:,jtype), zbias(:,:,:) )
219
220         jt=1
221         DO jobs = 1, sstdata%nsurf
222            IF ( sstdata%ntyp(jobs) == ibiastypes(jtype) ) THEN
223
224               zlam = sstdata%rlam(jobs)
225               zphi = sstdata%rphi(jobs)
226
227               CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,             &
228                  &                   zglam_tmp(:,:,jt), zgphi_tmp(:,:,jt), &
229                  &                   zmask_tmp(:,:,jt), zweig, zobsmask )
230
231               CALL obs_int_h2d( 1, 1, zweig, zbias(:,:,jt),  zext )
232
233               ! adjust sst with bias field
234               sstdata%robs(jobs,1) = sstdata%robs(jobs,1) - zext(1)
235
236               jt=jt+1
237
238            ENDIF
239         END DO
240
241         IF ( lwp ) THEN
242            WRITE(numout,*) 'Applied bias to obs of type: ',ibiastypes(jtype)
243            WRITE(numout,*) 'Max/min of bias: ', maxval(zbias(:,:,:)), minval(zbias(:,:,:))
244         ENDIF
245
246         !Deallocate arrays
247         DEALLOCATE( &
248         & igrdi_tmp, &
249         & igrdj_tmp, &
250         & zglam_tmp, &
251         & zgphi_tmp, &
252         & zmask_tmp, &
253         & zbias )     
254     
255      END DO !jtype
256
257      DEALLOCATE( &
258         & igrdi, &
259         & igrdj, &
260         & zglam, &
261         & zgphi, &
262         & zmask )
263
264      IF(lwp) THEN
265         WRITE(numout,*) " "
266         WRITE(numout,*) "SST bias correction applied successfully"
267         WRITE(numout,*) "Obs types: ",ibiastypes(:), &
268                              " have all been bias corrected\n"
269      ENDIF
270   END SUBROUTINE obs_app_sstbias
271 
272END MODULE obs_sstbias
Note: See TracBrowser for help on using the repository browser.