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_bias.F90 in NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_generic_obs/src/OCE/OBS/obs_bias.F90 @ 15187

Last change on this file since 15187 was 15187, checked in by dford, 11 months ago

Update treatment of SLA and POTM additional/extra variables.

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