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

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90 @ 2636

Last change on this file since 2636 was 2633, checked in by trackstand2, 13 years ago

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