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_readmdt.F90 in branches/dev_1784_OBS/NEMO/OPA_SRC/OBS – NEMO

source: branches/dev_1784_OBS/NEMO/OPA_SRC/OBS/obs_readmdt.F90 @ 2001

Last change on this file since 2001 was 2001, checked in by djlea, 14 years ago

Adding observation operator code
File size: 11.0 KB
Line 
1MODULE obs_readmdt
2   !!======================================================================
3   !!                       ***  MODULE obs_readmdt  ***
4   !! Observation diagnostics: Read the MDT for SLA data (skeleton for now)
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   obs_rea_mdt : Driver for reading MDT
9   !!----------------------------------------------------------------------
10
11   !! * Modules used   
12   USE par_kind, ONLY : &       ! Precision variables
13      & wp, &
14      & dp, &
15      & sp
16   USE par_oce, ONLY : &        ! Domain parameters
17      & jpi, &
18      & jpj, &
19      & jpim1
20   USE in_out_manager, ONLY : & ! I/O manager
21      & lwp,    &
22      & numout 
23   USE obs_surf_def             ! Surface observation definitions
24   USE dom_oce, ONLY : &        ! Domain variables
25      & tmask, &
26      & tmask_i, &
27      & e1t,   &
28      & e2t,   &
29      & gphit, &
30      & glamt
31   USE obs_const, ONLY : &
32      & obfillflt              ! Fillvalue
33   USE oce, ONLY : &           ! Model variables
34      & sshn
35   USE obs_inter_sup           ! Interpolation support routines
36   USE obs_inter_h2d           ! 2D interpolation
37   USE obs_utils               ! Various observation tools
38   USE lib_mpp, only: &        ! MPP routines
39      & lk_mpp, &
40      & mpp_sum
41   USE iom_nf90   
42   USE netcdf                   ! NetCDF library
43
44   IMPLICIT NONE
45
46   !! * Routine accessibility
47   PRIVATE
48
49   INTEGER, PUBLIC :: &
50      & nmsshc = 1        ! MDT correction scheme
51   REAL(wp), PUBLIC :: &
52      & mdtcorr = 1.61    ! User specified MDT correction
53   REAL(wp), PUBLIC :: &
54      & mdtcutoff = 65.0  ! MDT cutoff for computed correction
55   PUBLIC obs_rea_mdt     ! Read the MDT
56   PUBLIC obs_offset_mdt  ! Remove the offset between the model MDT and the
57                          ! used one
58
59CONTAINS
60
61   SUBROUTINE obs_rea_mdt( kslano, sladata, k2dint )
62      !!---------------------------------------------------------------------
63      !!
64      !!                   *** ROUTINE obs_rea_mdt ***
65      !!
66      !! ** Purpose : Read from file the MDT data (skeleton)
67      !!
68      !! ** Method  :
69      !!
70      !! ** Action  :
71      !!
72      !! References :
73      !!
74      !! History : 
75      !!      ! :  2007-03 (K. Mogensen) Initial skeleton version
76      !!      ! :  2007-04 (E. Remy) migration and improvement from OPAVAR
77      !!----------------------------------------------------------------------
78      !! * Modules used
79      USE iom
80
81      !! * Arguments
82      INTEGER, INTENT(IN) :: &
83         & kslano          ! Number of SLA Products
84      TYPE(obs_surf), DIMENSION(kslano), INTENT(INOUT) :: &
85         & sladata       ! SLA data
86      INTEGER, INTENT(IN) :: &
87         & k2dint
88
89      !! * Local declarations
90
91      CHARACTER(LEN=12), PARAMETER :: &
92         & cpname = 'obs_rea_mdt'
93      CHARACTER(LEN=20), PARAMETER :: &
94         & mdtname = 'slaReferenceLevel.nc'
95
96      INTEGER :: &
97         & jslano,  &      ! Data set loop variable
98         & jobs,    &      ! Obs loop variable
99         & jpimdt,  &      ! Number of grid point in latitude for the MDT
100         & jpjmdt,  &      ! Number of grid point in longitude for the MDT
101         & iico,    &      ! Grid point indicies
102         & ijco
103      INTEGER :: & 
104         & i_nx_id,     &  ! Index to read the NetCDF file
105         & i_ny_id,     &  !
106         & i_file_id,   &  !
107         & i_var_id,    &
108         & i_stat
109
110      REAL(wp), DIMENSION(jpi,jpj) :: & 
111         & z_mdt           ! Array to store the MDT values
112      REAL(wp), DIMENSION(1) :: &
113         & zext, &
114         & zobsmask
115      REAL(wp), DIMENSION(2,2,1) :: &
116         & zweig
117      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
118         & zmask, &
119         & zmdtl, &
120         & zglam, &
121         & zgphi
122      REAL(wp) :: &
123         & zlam, &
124         & zphi, &
125         & zfill
126      REAL(sp) :: zinfill
127      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
128         & igrdi, &
129         & igrdj
130      INTEGER :: &
131         & nummdt
132
133      IF(lwp)WRITE(numout,*) 
134      IF(lwp)WRITE(numout,*) ' obs_rea_mdt : '
135      IF(lwp)WRITE(numout,*) ' ------------- '
136      IF(lwp)WRITE(numout,*) '   Read MDT for referencing altimeter', &
137         &                   '   anomalies'
138
139      ! Open the file
140     
141      CALL iom_open( mdtname, nummdt )
142     
143      ! Get the MDT data
144
145      CALL iom_get( nummdt, jpdom_data, 'sossheig', z_mdt(:,:), 1 )
146
147      ! Close the file
148
149      CALL iom_close(nummdt)     
150     
151      ! Read in the fill value
152      zinfill = 0.0
153      i_stat = nf90_open( mdtname, nf90_nowrite, nummdt )
154      i_stat = nf90_inq_varid( nummdt, 'sossheig', i_var_id )
155      i_stat = nf90_get_att( nummdt, i_var_id, "_FillValue",zinfill)
156      zfill = zinfill
157      i_stat = nf90_close( nummdt )
158
159      ! Remove the offset between the MDT used with the sla and the model MDT
160
161      IF ( nmsshc == 1 .OR. nmsshc == 2 ) CALL obs_offset_mdt( z_mdt, zfill )
162
163      ! Intepolate the MDT already on the model grid at the observation point
164 
165      DO jslano = 1, kslano
166
167         ALLOCATE( &
168            & igrdi(2,2,sladata(jslano)%nsurf), &
169            & igrdj(2,2,sladata(jslano)%nsurf), &
170            & zglam(2,2,sladata(jslano)%nsurf), &
171            & zgphi(2,2,sladata(jslano)%nsurf), &
172            & zmask(2,2,sladata(jslano)%nsurf), &
173            & zmdtl(2,2,sladata(jslano)%nsurf)  &
174            & )
175         
176         DO jobs = 1, sladata(jslano)%nsurf
177
178            igrdi(1,1,jobs) = sladata(jslano)%mi(jobs)-1
179            igrdj(1,1,jobs) = sladata(jslano)%mj(jobs)-1
180            igrdi(1,2,jobs) = sladata(jslano)%mi(jobs)-1
181            igrdj(1,2,jobs) = sladata(jslano)%mj(jobs)
182            igrdi(2,1,jobs) = sladata(jslano)%mi(jobs)
183            igrdj(2,1,jobs) = sladata(jslano)%mj(jobs)-1
184            igrdi(2,2,jobs) = sladata(jslano)%mi(jobs)
185            igrdj(2,2,jobs) = sladata(jslano)%mj(jobs)
186
187         END DO
188
189         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
190            &                  igrdi, igrdj, glamt, zglam )
191         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
192            &                  igrdi, igrdj, gphit, zgphi )
193         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
194            &                  igrdi, igrdj, tmask(:,:,1), zmask )
195         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
196            &                  igrdi, igrdj, z_mdt, zmdtl )
197
198         DO jobs = 1, sladata(jslano)%nsurf
199           
200            zlam = sladata(jslano)%rlam(jobs)
201            zphi = sladata(jslano)%rphi(jobs)
202            iico = sladata(jslano)%mi(jobs)
203            ijco = sladata(jslano)%mj(jobs)
204           
205            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
206               &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
207               &                   zmask(:,:,jobs), zweig, zobsmask )
208           
209
210            IF ( z_mdt(iico-1,ijco-1) == zfill .OR. &
211                 z_mdt(iico-1,ijco  ) == zfill .OR. &
212                 z_mdt(iico  ,ijco-1) == zfill .OR. &
213                 z_mdt(iico  ,ijco  ) == zfill ) THEN
214               
215          sladata(jslano)%rext(jobs,2) = obfillflt
216               sladata(jslano)%nqc(jobs)    = 11                ! set qc flag for data with no mdt   
217       ELSE
218
219                CALL obs_int_h2d( 1, 1,      &
220                   &              zweig, zmdtl(:,:,jobs),  zext )
221           
222                sladata(jslano)%rext(jobs,2) = zext(1)
223
224            ENDIF
225
226         END DO
227         
228         DEALLOCATE( &
229            & igrdi, &
230            & igrdj, &
231            & zglam, &
232            & zgphi, &
233            & zmask, &
234            & zmdtl  &
235            & )
236
237      END DO
238
239   END SUBROUTINE obs_rea_mdt
240
241   SUBROUTINE obs_offset_mdt( mdt, zfill )
242
243      !!---------------------------------------------------------------------
244      !!
245      !!                   *** ROUTINE obs_offset_mdt ***
246      !!
247      !! ** Purpose : Compute a correction term for the MDT on the model grid
248      !!             !!!!! IF it is on the model grid
249      !!
250      !! ** Method  : Compute the mean difference between the model and the
251      !!              used MDT and remove the offset.
252      !!
253      !! ** Action  :
254      !!
255      !! References :
256      !!
257      !! History : 
258      !!      ! :  2007-04 (E. Remy) migration from OPAVAR
259      !!----------------------------------------------------------------------
260      !! * Modules used
261
262      !! * Arguments
263      REAL(wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
264         & mdt           ! MDT used on the model grid
265      REAL(wp), INTENT(IN) :: zfill 
266
267      !! * Local declarations
268      REAL(wp) :: & 
269         & zdxdy,       & 
270         & zarea,       & 
271         & zeta1,       & 
272         & zeta2,       & 
273         & zcorr_mdt,   & 
274         & zcorr_bcketa,&
275         & zcorr
276      REAL(wp), DIMENSION(jpi,jpj) :: &
277         & zpromsk
278      INTEGER :: &
279         & jj, &
280         & ji
281      CHARACTER(LEN=14), PARAMETER :: &
282         & cpname = 'obs_offset_mdt'
283   
284      !  Initialize the local mask, for domain projection
285      !  Also exclude mdt points which are set to missing data
286
287      DO ji = 1, jpi
288        DO jj = 1, jpj
289           zpromsk(ji,jj) = tmask_i(ji,jj)
290           IF (    ( gphit(ji,jj) .GT.  mdtcutoff ) &
291              &.OR.( gphit(ji,jj) .LT. -mdtcutoff ) &
292              &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
293              &        zpromsk(ji,jj) = 0.0
294        END DO
295      END DO 
296
297      ! Compute MSSH mean over [0,360] x [-mdtcutoff,mdtcutoff]
298
299      zarea = 0.0
300      zeta1 = 0.0
301      zeta2 = 0.0
302
303      DO jj = 1, jpj
304         DO ji = 1, jpi
305          zdxdy = e1t(ji,jj) * e2t(ji,jj) * zpromsk(ji,jj)
306          zarea = zarea + zdxdy
307          zeta1 = zeta1 + mdt(ji,jj) * zdxdy
308          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
309        END DO     
310      END DO
311
312      IF( lk_mpp) CALL mpp_sum( zeta1 )
313      IF( lk_mpp) CALL mpp_sum( zeta2 )
314      IF( lk_mpp) CALL mpp_sum( zarea )
315     
316      zcorr_mdt = zeta1 / zarea
317      zcorr_bcketa  = zeta2 / zarea
318
319      !  Define correction term
320
321      zcorr = zcorr_mdt - zcorr_bcketa
322
323      !  Correct spatial mean of the MSSH
324
325      IF ( nmsshc == 1 ) mdt(:,:) = mdt(:,:) - zcorr 
326
327      ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT
328
329      IF ( nmsshc == 2 ) mdt(:,:) = mdt(:,:) - mdtcorr
330
331      IF(lwp) THEN
332         WRITE(numout,*)
333         WRITE(numout,*) ' obs_readmdt : mdtcutoff     = ', mdtcutoff
334         WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
335         WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
336         WRITE(numout,*) '               zcorr         = ', zcorr
337         WRITE(numout,*) '               nmsshc        = ', nmsshc
338         WRITE(numout,*) 
339      ENDIF
340
341      IF ( nmsshc == 0 ) WRITE(numout,*) &
342         &                 '           MSSH correction is not applied'
343      IF ( nmsshc == 1 ) WRITE(numout,*) &
344         &                 '           MSSH correction is applied'
345      IF ( nmsshc == 2 ) WRITE(numout,*) &
346         &                 '           User defined MSSH correction' 
347
348
349   END SUBROUTINE obs_offset_mdt
350 
351END MODULE obs_readmdt
Note: See TracBrowser for help on using the repository browser.