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 @ 2638

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

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