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

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

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