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.
obs_sstbias.F90 in branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/OBS/obs_sstbias.F90 @ 4768

Last change on this file since 4768 was 4768, checked in by jwhile, 10 years ago

Adding code for bias correction to SST obs

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