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

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

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