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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
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   !! History :      ! 2007-03 (K. Mogensen) Initial skeleton version
7   !!                ! 2007-04 (E. Remy) migration and improvement from OPAVAR
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   obs_rea_mdt    : Driver for reading MDT
12   !!   obs_offset_mdt : Remove the offset between the model MDT and the used one
13   !!----------------------------------------------------------------------
14   USE par_kind         ! Precision variables
15   USE par_oce          ! Domain parameters
16   USE in_out_manager   ! I/O manager
17   USE obs_surf_def     ! Surface observation definitions
18   USE obs_inter_sup    ! Interpolation support routines
19   USE obs_inter_h2d    ! 2D interpolation
20   USE obs_utils        ! Various observation tools
21   USE iom_nf90         ! IOM NetCDF
22   USE netcdf           ! NetCDF library
23   USE lib_mpp          ! MPP library
24   USE dom_oce, ONLY : &                  ! Domain variables
25      &                    tmask, tmask_i, e1t, e2t, gphit, glamt
26   USE obs_const, ONLY :   obfillflt      ! Fillvalue
27   USE oce      , ONLY :   sshn           ! Model variables
28
29   IMPLICIT NONE
30   PRIVATE
31   
32   PUBLIC   obs_rea_mdt     ! called by ?
33   PUBLIC   obs_offset_mdt  ! called by ?
34
35   INTEGER , PUBLIC ::   nmsshc    = 1         ! MDT correction scheme
36   REAL(wp), PUBLIC ::   mdtcorr   = 1.61_wp   ! User specified MDT correction
37   REAL(wp), PUBLIC ::   mdtcutoff = 65.0_wp   ! MDT cutoff for computed correction
38
39   !! * Control permutation of array indices
40#  include "dom_oce_ftrans.h90"
41
42   !!----------------------------------------------------------------------
43   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
44   !! $Id$
45   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE obs_rea_mdt( kslano, sladata, k2dint )
50      !!---------------------------------------------------------------------
51      !!
52      !!                   *** ROUTINE obs_rea_mdt ***
53      !!
54      !! ** Purpose : Read from file the MDT data (skeleton)
55      !!
56      !! ** Method  :
57      !!
58      !! ** Action  :
59      !!----------------------------------------------------------------------
60      USE iom
61      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
62      USE wrk_nemo, ONLY:   z_mdt   => wrk_2d_1   ! Array to store the MDT values
63      USE wrk_nemo, ONLY:   mdtmask => wrk_2d_2   ! Array to store the mask for the MDT
64      !
65      INTEGER                          , INTENT(IN)    ::   kslano    ! Number of SLA Products
66      TYPE(obs_surf), DIMENSION(kslano), INTENT(inout) ::   sladata   ! SLA data
67      INTEGER                          , INTENT(in)    ::   k2dint    ! ?
68      !
69      CHARACTER(LEN=12), PARAMETER ::   cpname  = 'obs_rea_mdt'
70      CHARACTER(LEN=20), PARAMETER ::   mdtname = 'slaReferenceLevel.nc'
71
72      INTEGER ::   jslano              ! Data set loop variable
73      INTEGER ::   jobs                ! Obs loop variable
74      INTEGER ::   jpimdt, jpjmdt      ! Number of grid point in lat/lon for the MDT
75      INTEGER ::   iico, ijco          ! Grid point indicies
76      INTEGER ::   i_nx_id, i_ny_id, i_file_id, i_var_id, i_stat
77      INTEGER ::   nummdt
78      !
79      REAL(wp), DIMENSION(1)     ::   zext, zobsmask
80      REAL(wp), DIMENSION(2,2,1) ::   zweig
81      !
82      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zmask, zmdtl, zglam, zgphi
83      INTEGER , DIMENSION(:,:,:), ALLOCATABLE ::   igrdi, igrdj
84         
85      REAL(wp) :: zlam, zphi, zfill, zinfill    ! local scalar
86      !!----------------------------------------------------------------------
87
88      IF( wrk_in_use(2, 1,2) ) THEN
89         CALL ctl_stop('obs_rea_mdt : requested workspace array unavailable')   ;   RETURN
90      ENDIF
91
92      IF(lwp)WRITE(numout,*) 
93      IF(lwp)WRITE(numout,*) ' obs_rea_mdt : Read MDT for referencing altimeter anomalies'
94      IF(lwp)WRITE(numout,*) ' ------------- '
95
96      CALL iom_open( mdtname, nummdt )       ! Open the file
97      !                                      ! Get the MDT data
98      CALL iom_get ( nummdt, jpdom_data, 'sossheig', z_mdt(:,:), 1 )
99      CALL iom_close(nummdt)                 ! Close the file
100     
101      ! Read in the fill value
102      zinfill = 0.0
103      i_stat = nf90_open( mdtname, nf90_nowrite, nummdt )
104      i_stat = nf90_inq_varid( nummdt, 'sossheig', i_var_id )
105      i_stat = nf90_get_att( nummdt, i_var_id, "_FillValue",zinfill)
106      zfill = zinfill
107      i_stat = nf90_close( nummdt )
108
109      ! setup mask based on tmask and MDT mask
110      ! set mask to 0 where the MDT is set to fillvalue
111#if defined key_z_first
112      WHERE(z_mdt(:,:) /= zfill)   ;   mdtmask(:,:) = tmask_1(:,:)
113#else
114      WHERE(z_mdt(:,:) /= zfill)   ;   mdtmask(:,:) = tmask(:,:,1)
115#endif
116      ELSE WHERE                   ;   mdtmask(:,:) = 0
117      END WHERE
118
119      ! Remove the offset between the MDT used with the sla and the model MDT
120      IF( nmsshc == 1 .OR. nmsshc == 2 )   CALL obs_offset_mdt( z_mdt, zfill )
121
122      ! Intepolate the MDT already on the model grid at the observation point
123 
124      DO jslano = 1, kslano
125         ALLOCATE( &
126            & igrdi(2,2,sladata(jslano)%nsurf), &
127            & igrdj(2,2,sladata(jslano)%nsurf), &
128            & zglam(2,2,sladata(jslano)%nsurf), &
129            & zgphi(2,2,sladata(jslano)%nsurf), &
130            & zmask(2,2,sladata(jslano)%nsurf), &
131            & zmdtl(2,2,sladata(jslano)%nsurf)  &
132            & )
133         
134         DO jobs = 1, sladata(jslano)%nsurf
135
136            igrdi(1,1,jobs) = sladata(jslano)%mi(jobs)-1
137            igrdj(1,1,jobs) = sladata(jslano)%mj(jobs)-1
138            igrdi(1,2,jobs) = sladata(jslano)%mi(jobs)-1
139            igrdj(1,2,jobs) = sladata(jslano)%mj(jobs)
140            igrdi(2,1,jobs) = sladata(jslano)%mi(jobs)
141            igrdj(2,1,jobs) = sladata(jslano)%mj(jobs)-1
142            igrdi(2,2,jobs) = sladata(jslano)%mi(jobs)
143            igrdj(2,2,jobs) = sladata(jslano)%mj(jobs)
144
145         END DO
146
147         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, glamt  , zglam )
148         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, gphit  , zgphi )
149         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, mdtmask, zmask )
150         CALL obs_int_comm_2d( 2, 2, sladata(jslano)%nsurf, igrdi, igrdj, z_mdt  , zmdtl )
151
152         DO jobs = 1, sladata(jslano)%nsurf
153           
154            zlam = sladata(jslano)%rlam(jobs)
155            zphi = sladata(jslano)%rphi(jobs)
156
157            CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
158               &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
159               &                   zmask(:,:,jobs), zweig, zobsmask )
160           
161            CALL obs_int_h2d( 1, 1, zweig, zmdtl(:,:,jobs),  zext )
162 
163            sladata(jslano)%rext(jobs,2) = zext(1)
164
165! mark any masked data with a QC flag
166            IF( zobsmask(1) == 0 )   sladata(jslano)%nqc(jobs) = 11
167
168         END DO
169         
170         DEALLOCATE( &
171            & igrdi, &
172            & igrdj, &
173            & zglam, &
174            & zgphi, &
175            & zmask, &
176            & zmdtl  &
177            & )
178
179      END DO
180
181      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('obs_rea_mdt: failed to release workspace arrays')
182      !
183   END SUBROUTINE obs_rea_mdt
184
185
186   SUBROUTINE obs_offset_mdt( mdt, zfill )
187      !!---------------------------------------------------------------------
188      !!
189      !!                   *** ROUTINE obs_offset_mdt ***
190      !!
191      !! ** Purpose : Compute a correction term for the MDT on the model grid
192      !!             !!!!! IF it is on the model grid
193      !!
194      !! ** Method  : Compute the mean difference between the model and the
195      !!              used MDT and remove the offset.
196      !!
197      !! ** Action  :
198      !!----------------------------------------------------------------------
199      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
200      USE wrk_nemo, ONLY:   zpromsk => wrk_2d_3
201      !
202      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   mdt     ! MDT used on the model grid
203      REAL(wp)                    , INTENT(in   ) ::   zfill 
204      !
205      INTEGER  :: ji, jj
206      REAL(wp) :: zdxdy, zarea, zeta1, zeta2, zcorr_mdt, zcorr_bcketa, zcorr     ! local scalar
207      CHARACTER(LEN=14), PARAMETER ::   cpname = 'obs_offset_mdt'
208      !!----------------------------------------------------------------------
209
210      IF( wrk_in_use(2, 3) ) THEN
211         CALL ctl_stop('obs_offset_mdt: requested workspace array unavailable')   ;   RETURN
212      ENDIF
213
214      !  Initialize the local mask, for domain projection
215      !  Also exclude mdt points which are set to missing data
216
217      DO ji = 1, jpi
218        DO jj = 1, jpj
219           zpromsk(ji,jj) = tmask_i(ji,jj)
220           IF (    ( gphit(ji,jj) .GT.  mdtcutoff ) &
221              &.OR.( gphit(ji,jj) .LT. -mdtcutoff ) &
222              &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
223              &        zpromsk(ji,jj) = 0.0
224        END DO
225      END DO 
226
227      ! Compute MSSH mean over [0,360] x [-mdtcutoff,mdtcutoff]
228
229      zarea = 0.0
230      zeta1 = 0.0
231      zeta2 = 0.0
232
233      DO jj = 1, jpj
234         DO ji = 1, jpi
235          zdxdy = e1t(ji,jj) * e2t(ji,jj) * zpromsk(ji,jj)
236          zarea = zarea + zdxdy
237          zeta1 = zeta1 + mdt(ji,jj) * zdxdy
238          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
239        END DO     
240      END DO
241
242      IF( lk_mpp)   CALL mpp_sum( zeta1 )
243      IF( lk_mpp)   CALL mpp_sum( zeta2 )
244      IF( lk_mpp)   CALL mpp_sum( zarea )
245     
246      zcorr_mdt    = zeta1 / zarea
247      zcorr_bcketa = zeta2 / zarea
248
249      !  Define correction term
250
251      zcorr = zcorr_mdt - zcorr_bcketa
252
253      !  Correct spatial mean of the MSSH
254
255      IF( nmsshc == 1 )   mdt(:,:) = mdt(:,:) - zcorr 
256
257      ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT
258
259      IF( nmsshc == 2 )   mdt(:,:) = mdt(:,:) - mdtcorr
260
261      IF(lwp) THEN
262         WRITE(numout,*)
263         WRITE(numout,*) ' obs_readmdt : mdtcutoff     = ', mdtcutoff
264         WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
265         WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
266         WRITE(numout,*) '               zcorr         = ', zcorr
267         WRITE(numout,*) '               nmsshc        = ', nmsshc
268      ENDIF
269
270      IF ( nmsshc == 0 ) WRITE(numout,*) '           MSSH correction is not applied'
271      IF ( nmsshc == 1 ) WRITE(numout,*) '           MSSH correction is applied'
272      IF ( nmsshc == 2 ) WRITE(numout,*) '           User defined MSSH correction' 
273
274      IF( wrk_not_released(2, 3) )   CALL ctl_stop('obs_offset_mdt: failed to release workspace array')
275      !
276   END SUBROUTINE obs_offset_mdt
277 
278   !!======================================================================
279END MODULE obs_readmdt
Note: See TracBrowser for help on using the repository browser.