source: NEMO/trunk/src/OCE/OBS/obs_sstbias.F90

Last change on this file was 13286, checked in by smasson, 5 months ago

trunk: merge extra halos branch in trunk, see #2366

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