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/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update_addclim/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90 @ 11449

Last change on this file since 11449 was 11449, checked in by mattmartin, 5 years ago

Committed first version of changes to output climatology values at obs locations in the fdbk files.
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, e1t, e2t, 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 dia_obs_init
34   PUBLIC   obs_offset_mdt  ! called by obs_rea_mdt
35
36   INTEGER , PUBLIC :: nn_msshc    = 1         ! MDT correction scheme
37   REAL(wp), PUBLIC :: rn_mdtcorr   = 1.61_wp  ! User specified MDT correction
38   REAL(wp), PUBLIC :: rn_mdtcutoff = 65.0_wp  ! MDT cutoff for computed correction
39   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::  mdt_fld    !: Mean Dynamic Topography on the model grid
40   
41   !!----------------------------------------------------------------------
42   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
43   !! $Id$
44   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE obs_rea_mdt( sladata, k2dint )
49      !!---------------------------------------------------------------------
50      !!
51      !!                   *** ROUTINE obs_rea_mdt ***
52      !!
53      !! ** Purpose : Read from file the MDT data (skeleton)
54      !!
55      !! ** Method  :
56      !!
57      !! ** Action  :
58      !!----------------------------------------------------------------------
59      USE iom
60      !
61      TYPE(obs_surf), 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 ::   jobs                ! Obs loop variable
68      INTEGER ::   jpimdt, jpjmdt      ! Number of grid point in lat/lon for the MDT
69      INTEGER ::   iico, ijco          ! Grid point indicies
70      INTEGER ::   i_nx_id, i_ny_id, i_file_id, i_var_id, i_stat
71      INTEGER ::   nummdt
72      !
73      REAL(wp), DIMENSION(1)     ::   zext, zobsmask
74      REAL(wp), DIMENSION(2,2,1) ::   zweig
75      !
76      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zmask, zmdtl, zglam, zgphi
77      INTEGER , DIMENSION(:,:,:), ALLOCATABLE ::   igrdi, igrdj
78      !
79      REAL(wp), POINTER, DIMENSION(:,:) ::  z_mdt, mdtmask
80         
81      REAL(wp) :: zlam, zphi, zfill, zinfill    ! local scalar
82      !!----------------------------------------------------------------------
83
84      CALL wrk_alloc(jpi,jpj,z_mdt,mdtmask) 
85
86      IF(lwp)WRITE(numout,*) 
87      IF(lwp)WRITE(numout,*) ' obs_rea_mdt : Read MDT for referencing altimeter anomalies'
88      IF(lwp)WRITE(numout,*) ' ------------- '
89      CALL FLUSH(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      ! Set the MDT gridded field to be the one read in directly with no adjustments
97      ALLOCATE( mdt_fld(jpi,jpj) )
98      mdt_fld(:,:) = z_mdt(:,:)
99     
100      ! Read in the fill value
101      zinfill = 0.0
102      i_stat = nf90_open( mdtname, nf90_nowrite, nummdt )
103      i_stat = nf90_inq_varid( nummdt, 'sossheig', i_var_id )
104      i_stat = nf90_get_att( nummdt, i_var_id, "_FillValue",zinfill)
105      zfill = zinfill
106      i_stat = nf90_close( nummdt )
107
108      ! setup mask based on tmask and MDT mask
109      ! set mask to 0 where the MDT is set to fillvalue
110      WHERE(z_mdt(:,:) /= zfill)   ;   mdtmask(:,:) = tmask(:,:,1)
111      ELSE WHERE                   ;   mdtmask(:,:) = 0
112      END WHERE
113
114      ! Remove the offset between the MDT used with the sla and the model MDT
115      IF( nn_msshc == 1 .OR. nn_msshc == 2 ) &
116         & CALL obs_offset_mdt( jpi, jpj, z_mdt, zfill )
117
118      ! Intepolate the MDT already on the model grid at the observation point
119 
120      ALLOCATE( &
121         & igrdi(2,2,sladata%nsurf), &
122         & igrdj(2,2,sladata%nsurf), &
123         & zglam(2,2,sladata%nsurf), &
124         & zgphi(2,2,sladata%nsurf), &
125         & zmask(2,2,sladata%nsurf), &
126         & zmdtl(2,2,sladata%nsurf)  &
127         & )
128         
129      DO jobs = 1, sladata%nsurf
130
131         igrdi(1,1,jobs) = sladata%mi(jobs)-1
132         igrdj(1,1,jobs) = sladata%mj(jobs)-1
133         igrdi(1,2,jobs) = sladata%mi(jobs)-1
134         igrdj(1,2,jobs) = sladata%mj(jobs)
135         igrdi(2,1,jobs) = sladata%mi(jobs)
136         igrdj(2,1,jobs) = sladata%mj(jobs)-1
137         igrdi(2,2,jobs) = sladata%mi(jobs)
138         igrdj(2,2,jobs) = sladata%mj(jobs)
139
140      END DO
141
142      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, glamt  , zglam )
143      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, gphit  , zgphi )
144      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, mdtmask, zmask )
145      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, z_mdt  , zmdtl )
146
147      DO jobs = 1, sladata%nsurf
148           
149         zlam = sladata%rlam(jobs)
150         zphi = sladata%rphi(jobs)
151
152         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
153            &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
154            &                   zmask(:,:,jobs), zweig, zobsmask )
155           
156         CALL obs_int_h2d( 1, 1, zweig, zmdtl(:,:,jobs),  zext )
157 
158         sladata%rext(jobs,2) = zext(1)
159
160! mark any masked data with a QC flag
161         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15)
162
163         END DO
164         
165      DEALLOCATE( &
166         & igrdi, &
167         & igrdj, &
168         & zglam, &
169         & zgphi, &
170         & zmask, &
171         & zmdtl  &
172         & )
173
174      CALL wrk_dealloc(jpi,jpj,z_mdt,mdtmask) 
175      IF(lwp)WRITE(numout,*) ' ------------- '
176      !
177   END SUBROUTINE obs_rea_mdt
178
179
180   SUBROUTINE obs_offset_mdt( kpi, kpj, mdt, zfill )
181      !!---------------------------------------------------------------------
182      !!
183      !!                   *** ROUTINE obs_offset_mdt ***
184      !!
185      !! ** Purpose : Compute a correction term for the MDT on the model grid
186      !!             !!!!! IF it is on the model grid
187      !!
188      !! ** Method  : Compute the mean difference between the model and the
189      !!              used MDT and remove the offset.
190      !!
191      !! ** Action  :
192      !!----------------------------------------------------------------------
193      INTEGER, INTENT(IN) ::  kpi, kpj
194      REAL(wp), DIMENSION(kpi,kpj), INTENT(INOUT) ::   mdt     ! MDT used on the model grid
195      REAL(wp)                    , INTENT(IN   ) ::   zfill 
196      !
197      INTEGER  :: ji, jj
198      REAL(wp) :: zdxdy, zarea, zeta1, zeta2, zcorr_mdt, zcorr_bcketa, zcorr     ! local scalar
199      REAL(wp), POINTER, DIMENSION(:,:) :: zpromsk
200      CHARACTER(LEN=14), PARAMETER ::   cpname = 'obs_offset_mdt'
201      !!----------------------------------------------------------------------
202
203      CALL wrk_alloc( jpi,jpj, zpromsk )
204
205      !  Initialize the local mask, for domain projection
206      !  Also exclude mdt points which are set to missing data
207
208      DO ji = 1, jpi
209        DO jj = 1, jpj
210           zpromsk(ji,jj) = tmask_i(ji,jj)
211           IF (    ( gphit(ji,jj) .GT.  rn_mdtcutoff ) &
212              &.OR.( gphit(ji,jj) .LT. -rn_mdtcutoff ) &
213              &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
214              &        zpromsk(ji,jj) = 0.0
215        END DO
216      END DO 
217
218      ! Compute MSSH mean over [0,360] x [-rn_mdtcutoff,rn_mdtcutoff]
219
220      zarea = 0.0
221      zeta1 = 0.0
222      zeta2 = 0.0
223
224      DO jj = 1, jpj
225         DO ji = 1, jpi
226          zdxdy = e1t(ji,jj) * e2t(ji,jj) * zpromsk(ji,jj)
227          zarea = zarea + zdxdy
228          zeta1 = zeta1 + mdt(ji,jj) * zdxdy
229          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
230        END DO     
231      END DO
232
233      IF( lk_mpp)   CALL mpp_sum( zeta1 )
234      IF( lk_mpp)   CALL mpp_sum( zeta2 )
235      IF( lk_mpp)   CALL mpp_sum( zarea )
236     
237      zcorr_mdt    = zeta1 / zarea
238      zcorr_bcketa = zeta2 / zarea
239
240      !  Define correction term
241
242      zcorr = zcorr_mdt - zcorr_bcketa
243
244      !  Correct spatial mean of the MSSH
245
246      IF( nn_msshc == 1 )   mdt(:,:) = mdt(:,:) - zcorr 
247
248      ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT
249
250      IF( nn_msshc == 2 )   mdt(:,:) = mdt(:,:) - rn_mdtcorr
251
252      IF(lwp) THEN
253         WRITE(numout,*)
254         WRITE(numout,*) ' obs_readmdt : rn_mdtcutoff     = ', rn_mdtcutoff
255         WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
256         WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
257         WRITE(numout,*) '               zcorr         = ', zcorr
258         WRITE(numout,*) '               nn_msshc        = ', nn_msshc
259      ENDIF
260
261      IF ( nn_msshc == 0 ) WRITE(numout,*) '           MSSH correction is not applied'
262      IF ( nn_msshc == 1 ) WRITE(numout,*) '           MSSH correction is applied'
263      IF ( nn_msshc == 2 ) WRITE(numout,*) '           User defined MSSH correction' 
264
265      CALL wrk_dealloc( jpi,jpj, zpromsk )
266      !
267   END SUBROUTINE obs_offset_mdt
268 
269   !!======================================================================
270END MODULE obs_readmdt
Note: See TracBrowser for help on using the repository browser.