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/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90 @ 2287

Last change on this file since 2287 was 2287, checked in by smasson, 14 years ago

update licence of all NEMO files...
  • Property svn:keywords set to Id
File size: 11.1 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 :: nmsshc = 1        ! MDT correction scheme
50   REAL(wp), PUBLIC :: mdtcorr = 1.61   ! User specified MDT correction
51   REAL(wp), PUBLIC :: mdtcutoff = 65.0  ! MDT cutoff for computed correction
52   PUBLIC obs_rea_mdt     ! Read the MDT
53   PUBLIC obs_offset_mdt  ! Remove the offset between the model MDT and the
54                          ! used one
55
56   !!----------------------------------------------------------------------
57   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
58   !! $Id$
59   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
60   !!----------------------------------------------------------------------
61
62CONTAINS
63
64   SUBROUTINE obs_rea_mdt( kslano, sladata, k2dint )
65      !!---------------------------------------------------------------------
66      !!
67      !!                   *** ROUTINE obs_rea_mdt ***
68      !!
69      !! ** Purpose : Read from file the MDT data (skeleton)
70      !!
71      !! ** Method  :
72      !!
73      !! ** Action  :
74      !!
75      !! References :
76      !!
77      !! History : 
78      !!      ! :  2007-03 (K. Mogensen) Initial skeleton version
79      !!      ! :  2007-04 (E. Remy) migration and improvement from OPAVAR
80      !!----------------------------------------------------------------------
81      !! * Modules used
82      USE iom
83
84      !! * Arguments
85      INTEGER, INTENT(IN) :: kslano          ! Number of SLA Products
86      TYPE(obs_surf), DIMENSION(kslano), INTENT(INOUT) :: &
87         & sladata       ! SLA data
88      INTEGER, INTENT(IN) :: k2dint
89
90      !! * Local declarations
91
92      CHARACTER(LEN=12), PARAMETER :: &
93         & cpname = 'obs_rea_mdt'
94      CHARACTER(LEN=20), PARAMETER :: &
95         & mdtname = 'slaReferenceLevel.nc'
96
97      INTEGER :: jslano      ! Data set loop variable
98      INTEGER :: jobs        ! Obs loop variable
99      INTEGER :: jpimdt      ! Number of grid point in latitude for the MDT
100      INTEGER :: jpjmdt      ! Number of grid point in longitude for the MDT
101      INTEGER :: iico        ! Grid point indicies
102      INTEGER :: ijco 
103      INTEGER :: i_nx_id     ! Index to read the NetCDF file
104      INTEGER :: i_ny_id     !
105      INTEGER :: i_file_id   !
106      INTEGER :: i_var_id
107      INTEGER :: i_stat
108
109      REAL(wp), DIMENSION(jpi,jpj) :: & 
110         & z_mdt,       &  ! Array to store the MDT values
111         & mdtmask         ! Array to store the mask for the MDT
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         
123      REAL(wp) :: zlam
124      REAL(wp) :: zphi
125      REAL(wp) :: zfill
126      REAL(sp) :: zinfill
127      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: &
128         & igrdi, &
129         & igrdj
130      INTEGER :: nummdt
131
132      IF(lwp)WRITE(numout,*) 
133      IF(lwp)WRITE(numout,*) ' obs_rea_mdt : '
134      IF(lwp)WRITE(numout,*) ' ------------- '
135      IF(lwp)WRITE(numout,*) '   Read MDT for referencing altimeter', &
136         &                   '   anomalies'
137
138      ! Open the file
139     
140      CALL iom_open( mdtname, nummdt )
141     
142      ! Get the MDT data
143
144      CALL iom_get( nummdt, jpdom_data, 'sossheig', z_mdt(:,:), 1 )
145
146      ! Close the file
147
148      CALL iom_close(nummdt)     
149     
150      ! Read in the fill value
151      zinfill = 0.0
152      i_stat = nf90_open( mdtname, nf90_nowrite, nummdt )
153      i_stat = nf90_inq_varid( nummdt, 'sossheig', i_var_id )
154      i_stat = nf90_get_att( nummdt, i_var_id, "_FillValue",zinfill)
155      zfill = zinfill
156      i_stat = nf90_close( nummdt )
157
158! setup mask based on tmask and MDT mask
159! set mask to 0 where the MDT is set to fillvalue
160
161      WHERE(z_mdt(:,:) /= zfill)
162         mdtmask(:,:)=tmask(:,:,1)
163      ELSEWHERE
164         mdtmask(:,:)=0
165      END WHERE
166
167      ! Remove the offset between the MDT used with the sla and the model MDT
168
169      IF ( nmsshc == 1 .OR. nmsshc == 2 ) CALL obs_offset_mdt( z_mdt, zfill )
170
171      ! Intepolate the MDT already on the model grid at the observation point
172 
173      DO jslano = 1, kslano
174
175         ALLOCATE( &
176            & igrdi(2,2,sladata(jslano)%nsurf), &
177            & igrdj(2,2,sladata(jslano)%nsurf), &
178            & zglam(2,2,sladata(jslano)%nsurf), &
179            & zgphi(2,2,sladata(jslano)%nsurf), &
180            & zmask(2,2,sladata(jslano)%nsurf), &
181            & zmdtl(2,2,sladata(jslano)%nsurf)  &
182            & )
183         
184         DO jobs = 1, sladata(jslano)%nsurf
185
186            igrdi(1,1,jobs) = sladata(jslano)%mi(jobs)-1
187            igrdj(1,1,jobs) = sladata(jslano)%mj(jobs)-1
188            igrdi(1,2,jobs) = sladata(jslano)%mi(jobs)-1
189            igrdj(1,2,jobs) = sladata(jslano)%mj(jobs)
190            igrdi(2,1,jobs) = sladata(jslano)%mi(jobs)
191            igrdj(2,1,jobs) = sladata(jslano)%mj(jobs)-1
192            igrdi(2,2,jobs) = sladata(jslano)%mi(jobs)
193            igrdj(2,2,jobs) = sladata(jslano)%mj(jobs)
194
195         END DO
196
197         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
198            &                  igrdi, igrdj, glamt, zglam )
199         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
200            &                  igrdi, igrdj, gphit, zgphi )
201         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
202            &                  igrdi, igrdj, mdtmask, zmask )
203         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, &
204            &                  igrdi, igrdj, z_mdt, zmdtl )
205
206         DO jobs = 1, sladata(jslano)%nsurf
207           
208            zlam = sladata(jslano)%rlam(jobs)
209            zphi = sladata(jslano)%rphi(jobs)
210
211            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
212               &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
213               &                   zmask(:,:,jobs), zweig, zobsmask )
214           
215            CALL obs_int_h2d( 1, 1,      &
216                   &              zweig, zmdtl(:,:,jobs),  zext )
217 
218            sladata(jslano)%rext(jobs,2) = zext(1)
219
220! mark any masked data with a QC flag
221            IF ( zobsmask(1) == 0 ) sladata(jslano)%nqc(jobs) = 11
222
223         END DO
224         
225         DEALLOCATE( &
226            & igrdi, &
227            & igrdj, &
228            & zglam, &
229            & zgphi, &
230            & zmask, &
231            & zmdtl  &
232            & )
233
234      END DO
235
236   END SUBROUTINE obs_rea_mdt
237
238   SUBROUTINE obs_offset_mdt( mdt, zfill )
239
240      !!---------------------------------------------------------------------
241      !!
242      !!                   *** ROUTINE obs_offset_mdt ***
243      !!
244      !! ** Purpose : Compute a correction term for the MDT on the model grid
245      !!             !!!!! IF it is on the model grid
246      !!
247      !! ** Method  : Compute the mean difference between the model and the
248      !!              used MDT and remove the offset.
249      !!
250      !! ** Action  :
251      !!
252      !! References :
253      !!
254      !! History : 
255      !!      ! :  2007-04 (E. Remy) migration from OPAVAR
256      !!----------------------------------------------------------------------
257      !! * Modules used
258
259      !! * Arguments
260      REAL(wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
261         & mdt           ! MDT used on the model grid
262      REAL(wp), INTENT(IN) :: zfill 
263
264      !! * Local declarations
265      REAL(wp) :: zdxdy
266      REAL(wp) :: zarea
267      REAL(wp) :: zeta1
268      REAL(wp) :: zeta2
269      REAL(wp) :: zcorr_mdt 
270      REAL(wp) :: zcorr_bcketa
271      REAL(wp) :: zcorr
272      REAL(wp), DIMENSION(jpi,jpj) :: zpromsk
273      INTEGER :: jj
274      INTEGER :: ji
275      CHARACTER(LEN=14), PARAMETER :: &
276         & cpname = 'obs_offset_mdt'
277   
278      !  Initialize the local mask, for domain projection
279      !  Also exclude mdt points which are set to missing data
280
281      DO ji = 1, jpi
282        DO jj = 1, jpj
283           zpromsk(ji,jj) = tmask_i(ji,jj)
284           IF (    ( gphit(ji,jj) .GT.  mdtcutoff ) &
285              &.OR.( gphit(ji,jj) .LT. -mdtcutoff ) &
286              &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
287              &        zpromsk(ji,jj) = 0.0
288        END DO
289      END DO 
290
291      ! Compute MSSH mean over [0,360] x [-mdtcutoff,mdtcutoff]
292
293      zarea = 0.0
294      zeta1 = 0.0
295      zeta2 = 0.0
296
297      DO jj = 1, jpj
298         DO ji = 1, jpi
299          zdxdy = e1t(ji,jj) * e2t(ji,jj) * zpromsk(ji,jj)
300          zarea = zarea + zdxdy
301          zeta1 = zeta1 + mdt(ji,jj) * zdxdy
302          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
303        END DO     
304      END DO
305
306      IF( lk_mpp) CALL mpp_sum( zeta1 )
307      IF( lk_mpp) CALL mpp_sum( zeta2 )
308      IF( lk_mpp) CALL mpp_sum( zarea )
309     
310      zcorr_mdt = zeta1 / zarea
311      zcorr_bcketa  = zeta2 / zarea
312
313      !  Define correction term
314
315      zcorr = zcorr_mdt - zcorr_bcketa
316
317      !  Correct spatial mean of the MSSH
318
319      IF ( nmsshc == 1 ) mdt(:,:) = mdt(:,:) - zcorr 
320
321      ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT
322
323      IF ( nmsshc == 2 ) mdt(:,:) = mdt(:,:) - mdtcorr
324
325      IF(lwp) THEN
326         WRITE(numout,*)
327         WRITE(numout,*) ' obs_readmdt : mdtcutoff     = ', mdtcutoff
328         WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
329         WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
330         WRITE(numout,*) '               zcorr         = ', zcorr
331         WRITE(numout,*) '               nmsshc        = ', nmsshc
332         WRITE(numout,*) 
333      ENDIF
334
335      IF ( nmsshc == 0 ) WRITE(numout,*) &
336         &                 '           MSSH correction is not applied'
337      IF ( nmsshc == 1 ) WRITE(numout,*) &
338         &                 '           MSSH correction is applied'
339      IF ( nmsshc == 2 ) WRITE(numout,*) &
340         &                 '           User defined MSSH correction' 
341
342
343   END SUBROUTINE obs_offset_mdt
344 
345END MODULE obs_readmdt
Note: See TracBrowser for help on using the repository browser.