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

Last change on this file since 3849 was 3849, checked in by trackstand2, 11 years ago

Merge branch 'partitioner'
  • 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      ELSEWHERE                   ;   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.