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 NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/OBS/obs_readmdt.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago

Mixed precision version, tested up to 30 years on ORCA2.
  • Property svn:keywords set to Id
File size: 9.9 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, e1e2t, gphit, glamt
26   USE obs_const, ONLY :   obfillflt      ! Fillvalue
27   USE oce      , ONLY :   ssh            ! Model variables
28
29   IMPLICIT NONE
30   PRIVATE
31   
32   PUBLIC   obs_rea_mdt     ! called by dia_obs_init
33   PUBLIC   obs_offset_mdt  ! called by obs_rea_mdt
34
35   INTEGER , PUBLIC :: nn_msshc    = 1         ! MDT correction scheme
36   REAL(wp), PUBLIC :: rn_mdtcorr   = 1.61_wp  ! User specified MDT correction
37   REAL(wp), PUBLIC :: rn_mdtcutoff = 65.0_wp  ! MDT cutoff for computed correction
38
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41#  include "single_precision_substitute.h90"
42
43   !!----------------------------------------------------------------------
44   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
45   !! $Id$
46   !! Software governed by the CeCILL license (see ./LICENSE)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE obs_rea_mdt( sladata, k2dint, Kmm )
51      !!---------------------------------------------------------------------
52      !!
53      !!                   *** ROUTINE obs_rea_mdt ***
54      !!
55      !! ** Purpose : Read from file the MDT data (skeleton)
56      !!
57      !! ** Method  :
58      !!
59      !! ** Action  :
60      !!----------------------------------------------------------------------
61      USE iom
62      !
63      TYPE(obs_surf), INTENT(inout) ::   sladata   ! SLA data
64      INTEGER       , INTENT(in)    ::   k2dint    ! ?
65      INTEGER       , INTENT(in)    ::   Kmm       ! ?
66      !
67      CHARACTER(LEN=12), PARAMETER ::   cpname  = 'obs_rea_mdt'
68      CHARACTER(LEN=20), PARAMETER ::   mdtname = 'slaReferenceLevel.nc'
69
70      INTEGER ::   jobs                ! Obs loop variable
71      INTEGER ::   jpimdt, jpjmdt      ! Number of grid point in lat/lon for the MDT
72      INTEGER ::   iico, ijco          ! Grid point indicies
73      INTEGER ::   i_nx_id, i_ny_id, i_file_id, i_var_id, i_stat
74      INTEGER ::   nummdt
75      !
76      REAL(dp), DIMENSION(1)     ::   zext, zobsmask
77      REAL(dp), DIMENSION(2,2,1) ::   zweig
78      !
79      REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::   zmask, zmdtl, zglam, zgphi
80      INTEGER , DIMENSION(:,:,:), ALLOCATABLE ::   igrdi, igrdj
81      !
82      REAL(dp), DIMENSION(jpi,jpj) ::  z_mdt, mdtmask
83         
84      REAL(dp) :: zlam, zphi, zfill, zinfill    ! local scalar
85      !!----------------------------------------------------------------------
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      CALL FLUSH(numout)
91
92      CALL iom_open( mdtname, nummdt )       ! Open the file
93      !                                      ! Get the MDT data
94      CALL iom_get ( nummdt, jpdom_global, 'sossheig', z_mdt(:,:) )
95      CALL iom_close(nummdt)                 ! Close the file
96     
97      ! Read in the fill value
98      zinfill = 0.0
99      i_stat = nf90_open( mdtname, nf90_nowrite, nummdt )
100      i_stat = nf90_inq_varid( nummdt, 'sossheig', i_var_id )
101      i_stat = nf90_get_att( nummdt, i_var_id, "_FillValue",zinfill)
102      zfill = zinfill
103      i_stat = nf90_close( nummdt )
104
105      ! setup mask based on tmask and MDT mask
106      ! set mask to 0 where the MDT is set to fillvalue
107      WHERE(z_mdt(:,:) /= zfill)   ;   mdtmask(:,:) = tmask(:,:,1)
108      ELSE WHERE                   ;   mdtmask(:,:) = 0
109      END WHERE
110
111      ! Remove the offset between the MDT used with the sla and the model MDT
112      IF( nn_msshc == 1 .OR. nn_msshc == 2 ) &
113         & CALL obs_offset_mdt( jpi, jpj, z_mdt, zfill, Kmm )
114
115      ! Intepolate the MDT already on the model grid at the observation point
116 
117      ALLOCATE( &
118         & igrdi(2,2,sladata%nsurf), &
119         & igrdj(2,2,sladata%nsurf), &
120         & zglam(2,2,sladata%nsurf), &
121         & zgphi(2,2,sladata%nsurf), &
122         & zmask(2,2,sladata%nsurf), &
123         & zmdtl(2,2,sladata%nsurf)  &
124         & )
125         
126      DO jobs = 1, sladata%nsurf
127
128         igrdi(1,1,jobs) = sladata%mi(jobs)-1
129         igrdj(1,1,jobs) = sladata%mj(jobs)-1
130         igrdi(1,2,jobs) = sladata%mi(jobs)-1
131         igrdj(1,2,jobs) = sladata%mj(jobs)
132         igrdi(2,1,jobs) = sladata%mi(jobs)
133         igrdj(2,1,jobs) = sladata%mj(jobs)-1
134         igrdi(2,2,jobs) = sladata%mi(jobs)
135         igrdj(2,2,jobs) = sladata%mj(jobs)
136
137      END DO
138
139      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, glamt  , zglam )
140      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, gphit  , zgphi )
141      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, mdtmask, zmask )
142      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, z_mdt  , zmdtl )
143
144      DO jobs = 1, sladata%nsurf
145           
146         zlam = sladata%rlam(jobs)
147         zphi = sladata%rphi(jobs)
148
149         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
150            &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
151            &                   zmask(:,:,jobs), zweig, zobsmask )
152           
153         CALL obs_int_h2d( 1, 1, zweig, zmdtl(:,:,jobs),  zext )
154 
155         sladata%rext(jobs,2) = zext(1)
156
157! mark any masked data with a QC flag
158         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = IBSET(sladata%nqc(jobs),15)
159
160         END DO
161         
162      DEALLOCATE( &
163         & igrdi, &
164         & igrdj, &
165         & zglam, &
166         & zgphi, &
167         & zmask, &
168         & zmdtl  &
169         & )
170
171      IF(lwp)WRITE(numout,*) ' ------------- '
172      !
173   END SUBROUTINE obs_rea_mdt
174
175
176   SUBROUTINE obs_offset_mdt( kpi, kpj, mdt, zfill, Kmm )
177      !!---------------------------------------------------------------------
178      !!
179      !!                   *** ROUTINE obs_offset_mdt ***
180      !!
181      !! ** Purpose : Compute a correction term for the MDT on the model grid
182      !!             !!!!! IF it is on the model grid
183      !!
184      !! ** Method  : Compute the mean difference between the model and the
185      !!              used MDT and remove the offset.
186      !!
187      !! ** Action  :
188      !!----------------------------------------------------------------------
189      INTEGER, INTENT(IN) ::  kpi, kpj
190      INTEGER, INTENT(IN) ::  Kmm
191      REAL(dp), DIMENSION(kpi,kpj), INTENT(INOUT) ::   mdt     ! MDT used on the model grid
192      REAL(dp)                    , INTENT(IN   ) ::   zfill 
193      !
194      INTEGER  :: ji, jj
195      REAL(dp) :: zdxdy, zarea, zeta1, zeta2, zcorr_mdt, zcorr_bcketa, zcorr     ! local scalar
196      REAL(dp), DIMENSION(jpi,jpj) :: zpromsk
197      CHARACTER(LEN=14), PARAMETER ::   cpname = 'obs_offset_mdt'
198      !!----------------------------------------------------------------------
199
200      !  Initialize the local mask, for domain projection
201      !  Also exclude mdt points which are set to missing data
202
203      DO ji = 1, jpi
204        DO jj = 1, jpj
205           zpromsk(ji,jj) = tmask_i(ji,jj)
206           IF (    ( gphit(ji,jj) .GT.  rn_mdtcutoff ) &
207              &.OR.( gphit(ji,jj) .LT. -rn_mdtcutoff ) &
208              &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
209              &        zpromsk(ji,jj) = 0.0
210        END DO
211      END DO 
212
213      ! Compute MSSH mean over [0,360] x [-rn_mdtcutoff,rn_mdtcutoff]
214
215      zarea = 0.0
216      zeta1 = 0.0
217      zeta2 = 0.0
218
219      DO_2D( 1, 1, 1, 1 )
220       zdxdy = e1e2t(ji,jj) * zpromsk(ji,jj)
221       zarea = zarea + zdxdy
222       zeta1 = zeta1 + mdt(ji,jj) * zdxdy
223       zeta2 = zeta2 + ssh(ji,jj,Kmm) * zdxdy
224      END_2D
225
226      CALL mpp_sum( 'obs_readmdt', zeta1 )
227      CALL mpp_sum( 'obs_readmdt', zeta2 )
228      CALL mpp_sum( 'obs_readmdt', zarea )
229     
230      zcorr_mdt    = zeta1 / zarea
231      zcorr_bcketa = zeta2 / zarea
232
233      !  Define correction term
234
235      zcorr = zcorr_mdt - zcorr_bcketa
236
237      !  Correct spatial mean of the MSSH
238
239      IF( nn_msshc == 1 )   mdt(:,:) = mdt(:,:) - zcorr 
240
241      ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT
242
243      IF( nn_msshc == 2 )   mdt(:,:) = mdt(:,:) - rn_mdtcorr
244
245      IF(lwp) THEN
246         WRITE(numout,*)
247         WRITE(numout,*) ' obs_readmdt : rn_mdtcutoff     = ', rn_mdtcutoff
248         WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
249         WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
250         WRITE(numout,*) '               zcorr         = ', zcorr
251         WRITE(numout,*) '               nn_msshc        = ', nn_msshc
252      ENDIF
253
254      IF ( nn_msshc == 0 ) WRITE(numout,*) '           MSSH correction is not applied'
255      IF ( nn_msshc == 1 ) WRITE(numout,*) '           MSSH correction is applied'
256      IF ( nn_msshc == 2 ) WRITE(numout,*) '           User defined MSSH correction' 
257
258      !
259   END SUBROUTINE obs_offset_mdt
260 
261   !!======================================================================
262END MODULE obs_readmdt
Note: See TracBrowser for help on using the repository browser.