source: branches/UKMO/dev_r5785_SSS_obsoper/NEMOGCM/NEMO/OPA_SRC/OBS/obs_readmdt.F90 @ 7773

Last change on this file since 7773 was 7773, checked in by mattmartin, 4 years ago

Committing updates after doing the following:

  • merging the branch dev_r4650_general_vert_coord_obsoper@7763 into this branch
  • updating it so that the following OBS changes were implemented correctly on top of the simplification changes:
    • generalised vertical coordinate for profile obs. This was done so that is now the default option.
    • sst bias correction implemented with the new simplified obs code.
    • included the biogeochemical obs types int he new simplified obs code.
    • included the changes to exclude obs in the boundary for limited area models
    • included other changes for the efficiency of the obs operator to remove global arrays.
File size: 10.0 KB
RevLine 
[2128]1MODULE obs_readmdt
2   !!======================================================================
3   !!                       ***  MODULE obs_readmdt  ***
4   !! Observation diagnostics: Read the MDT for SLA data (skeleton for now)
5   !!======================================================================
[2715]6   !! History :      ! 2007-03 (K. Mogensen) Initial skeleton version
7   !!                ! 2007-04 (E. Remy) migration and improvement from OPAVAR
8   !!----------------------------------------------------------------------
[2128]9
10   !!----------------------------------------------------------------------
[2715]11   !!   obs_rea_mdt    : Driver for reading MDT
12   !!   obs_offset_mdt : Remove the offset between the model MDT and the used one
[2128]13   !!----------------------------------------------------------------------
[3294]14   USE wrk_nemo         ! Memory Allocation
[2715]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
[2128]29
30   IMPLICIT NONE
31   PRIVATE
[2715]32   
[5682]33   PUBLIC   obs_rea_mdt     ! called by dia_obs_init
34   PUBLIC   obs_offset_mdt  ! called by obs_rea_mdt
[2128]35
[5682]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
[2128]39
[2287]40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id$
[2715]43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[2287]44   !!----------------------------------------------------------------------
[2128]45CONTAINS
46
[5682]47   SUBROUTINE obs_rea_mdt( sladata, k2dint )
[2128]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
[2715]59      !
[5682]60      TYPE(obs_surf), INTENT(inout) ::   sladata   ! SLA data
61      INTEGER       , INTENT(in)    ::   k2dint    ! ?
[2715]62      !
63      CHARACTER(LEN=12), PARAMETER ::   cpname  = 'obs_rea_mdt'
64      CHARACTER(LEN=20), PARAMETER ::   mdtname = 'slaReferenceLevel.nc'
[2128]65
[2715]66      INTEGER ::   jobs                ! Obs loop variable
67      INTEGER ::   jpimdt, jpjmdt      ! Number of grid point in lat/lon for the MDT
68      INTEGER ::   iico, ijco          ! Grid point indicies
69      INTEGER ::   i_nx_id, i_ny_id, i_file_id, i_var_id, i_stat
70      INTEGER ::   nummdt
71      !
72      REAL(wp), DIMENSION(1)     ::   zext, zobsmask
73      REAL(wp), DIMENSION(2,2,1) ::   zweig
74      !
75      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   zmask, zmdtl, zglam, zgphi
76      INTEGER , DIMENSION(:,:,:), ALLOCATABLE ::   igrdi, igrdj
[3294]77      !
78      REAL(wp), POINTER, DIMENSION(:,:) ::  z_mdt, mdtmask
[2715]79         
80      REAL(wp) :: zlam, zphi, zfill, zinfill    ! local scalar
81      !!----------------------------------------------------------------------
[2128]82
[3294]83      CALL wrk_alloc(jpi,jpj,z_mdt,mdtmask) 
[2128]84
85      IF(lwp)WRITE(numout,*) 
[2715]86      IF(lwp)WRITE(numout,*) ' obs_rea_mdt : Read MDT for referencing altimeter anomalies'
[2128]87      IF(lwp)WRITE(numout,*) ' ------------- '
[5682]88      CALL FLUSH(numout)
[2128]89
[2715]90      CALL iom_open( mdtname, nummdt )       ! Open the file
91      !                                      ! Get the MDT data
92      CALL iom_get ( nummdt, jpdom_data, 'sossheig', z_mdt(:,:), 1 )
93      CALL iom_close(nummdt)                 ! Close the file
[2128]94     
95      ! Read in the fill value
96      zinfill = 0.0
97      i_stat = nf90_open( mdtname, nf90_nowrite, nummdt )
98      i_stat = nf90_inq_varid( nummdt, 'sossheig', i_var_id )
99      i_stat = nf90_get_att( nummdt, i_var_id, "_FillValue",zinfill)
100      zfill = zinfill
101      i_stat = nf90_close( nummdt )
102
[2715]103      ! setup mask based on tmask and MDT mask
104      ! set mask to 0 where the MDT is set to fillvalue
105      WHERE(z_mdt(:,:) /= zfill)   ;   mdtmask(:,:) = tmask(:,:,1)
106      ELSE WHERE                   ;   mdtmask(:,:) = 0
[2128]107      END WHERE
108
109      ! Remove the offset between the MDT used with the sla and the model MDT
[5704]110      IF( nn_msshc == 1 .OR. nn_msshc == 2 ) &
111         & CALL obs_offset_mdt( jpi, jpj, z_mdt, zfill )
[2128]112
113      ! Intepolate the MDT already on the model grid at the observation point
114 
[5682]115      ALLOCATE( &
116         & igrdi(2,2,sladata%nsurf), &
117         & igrdj(2,2,sladata%nsurf), &
118         & zglam(2,2,sladata%nsurf), &
119         & zgphi(2,2,sladata%nsurf), &
120         & zmask(2,2,sladata%nsurf), &
121         & zmdtl(2,2,sladata%nsurf)  &
122         & )
[2128]123         
[5682]124      DO jobs = 1, sladata%nsurf
[2128]125
[5682]126         igrdi(1,1,jobs) = sladata%mi(jobs)-1
127         igrdj(1,1,jobs) = sladata%mj(jobs)-1
128         igrdi(1,2,jobs) = sladata%mi(jobs)-1
129         igrdj(1,2,jobs) = sladata%mj(jobs)
130         igrdi(2,1,jobs) = sladata%mi(jobs)
131         igrdj(2,1,jobs) = sladata%mj(jobs)-1
132         igrdi(2,2,jobs) = sladata%mi(jobs)
133         igrdj(2,2,jobs) = sladata%mj(jobs)
[2128]134
[5682]135      END DO
[2128]136
[5704]137      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, glamt  , zglam )
138      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, gphit  , zgphi )
139      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, mdtmask, zmask )
140      CALL obs_int_comm_2d( 2, 2, sladata%nsurf, jpi, jpj, igrdi, igrdj, z_mdt  , zmdtl )
[2128]141
[5682]142      DO jobs = 1, sladata%nsurf
[2128]143           
[5682]144         zlam = sladata%rlam(jobs)
145         zphi = sladata%rphi(jobs)
[2128]146
[5682]147         CALL obs_int_h2d_init( 1, 1, k2dint, zlam, zphi,         &
148            &                   zglam(:,:,jobs), zgphi(:,:,jobs), &
149            &                   zmask(:,:,jobs), zweig, zobsmask )
[2128]150           
[5682]151         CALL obs_int_h2d( 1, 1, zweig, zmdtl(:,:,jobs),  zext )
[2128]152 
[5682]153         sladata%rext(jobs,2) = zext(1)
[2128]154
155! mark any masked data with a QC flag
[5682]156         IF( zobsmask(1) == 0 )   sladata%nqc(jobs) = 11
[2128]157
158         END DO
159         
[5682]160      DEALLOCATE( &
161         & igrdi, &
162         & igrdj, &
163         & zglam, &
164         & zgphi, &
165         & zmask, &
166         & zmdtl  &
167         & )
[2128]168
[3294]169      CALL wrk_dealloc(jpi,jpj,z_mdt,mdtmask) 
[5682]170      IF(lwp)WRITE(numout,*) ' ------------- '
[2715]171      !
[2128]172   END SUBROUTINE obs_rea_mdt
173
[2715]174
[5704]175   SUBROUTINE obs_offset_mdt( kpi, kpj, mdt, zfill )
[2128]176      !!---------------------------------------------------------------------
177      !!
178      !!                   *** ROUTINE obs_offset_mdt ***
179      !!
180      !! ** Purpose : Compute a correction term for the MDT on the model grid
181      !!             !!!!! IF it is on the model grid
182      !!
183      !! ** Method  : Compute the mean difference between the model and the
184      !!              used MDT and remove the offset.
185      !!
186      !! ** Action  :
187      !!----------------------------------------------------------------------
[5704]188      INTEGER, INTENT(IN) ::  kpi, kpj
189      REAL(wp), DIMENSION(kpi,kpj), INTENT(INOUT) ::   mdt     ! MDT used on the model grid
190      REAL(wp)                    , INTENT(IN   ) ::   zfill 
[2715]191      !
192      INTEGER  :: ji, jj
193      REAL(wp) :: zdxdy, zarea, zeta1, zeta2, zcorr_mdt, zcorr_bcketa, zcorr     ! local scalar
[3294]194      REAL(wp), POINTER, DIMENSION(:,:) :: zpromsk
[2715]195      CHARACTER(LEN=14), PARAMETER ::   cpname = 'obs_offset_mdt'
196      !!----------------------------------------------------------------------
[2128]197
[3294]198      CALL wrk_alloc( jpi,jpj, zpromsk )
[2128]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)
[5682]206           IF (    ( gphit(ji,jj) .GT.  rn_mdtcutoff ) &
207              &.OR.( gphit(ji,jj) .LT. -rn_mdtcutoff ) &
[2128]208              &.OR.( mdt(ji,jj) .EQ. zfill ) ) &
209              &        zpromsk(ji,jj) = 0.0
210        END DO
211      END DO 
212
[5682]213      ! Compute MSSH mean over [0,360] x [-rn_mdtcutoff,rn_mdtcutoff]
[2128]214
215      zarea = 0.0
216      zeta1 = 0.0
217      zeta2 = 0.0
218
219      DO jj = 1, jpj
220         DO ji = 1, jpi
221          zdxdy = e1t(ji,jj) * e2t(ji,jj) * zpromsk(ji,jj)
222          zarea = zarea + zdxdy
223          zeta1 = zeta1 + mdt(ji,jj) * zdxdy
224          zeta2 = zeta2 + sshn (ji,jj) * zdxdy
225        END DO     
226      END DO
227
[2715]228      IF( lk_mpp)   CALL mpp_sum( zeta1 )
229      IF( lk_mpp)   CALL mpp_sum( zeta2 )
230      IF( lk_mpp)   CALL mpp_sum( zarea )
[2128]231     
[2715]232      zcorr_mdt    = zeta1 / zarea
233      zcorr_bcketa = zeta2 / zarea
[2128]234
235      !  Define correction term
236
237      zcorr = zcorr_mdt - zcorr_bcketa
238
239      !  Correct spatial mean of the MSSH
240
[5682]241      IF( nn_msshc == 1 )   mdt(:,:) = mdt(:,:) - zcorr 
[2128]242
243      ! User defined value : 1.6 m for the Rio MDT compared to ORCA2 MDT
244
[5682]245      IF( nn_msshc == 2 )   mdt(:,:) = mdt(:,:) - rn_mdtcorr
[2128]246
247      IF(lwp) THEN
248         WRITE(numout,*)
[5682]249         WRITE(numout,*) ' obs_readmdt : rn_mdtcutoff     = ', rn_mdtcutoff
[2128]250         WRITE(numout,*) ' -----------   zcorr_mdt     = ', zcorr_mdt
251         WRITE(numout,*) '               zcorr_bcketa  = ', zcorr_bcketa
252         WRITE(numout,*) '               zcorr         = ', zcorr
[5682]253         WRITE(numout,*) '               nn_msshc        = ', nn_msshc
[2128]254      ENDIF
255
[5682]256      IF ( nn_msshc == 0 ) WRITE(numout,*) '           MSSH correction is not applied'
257      IF ( nn_msshc == 1 ) WRITE(numout,*) '           MSSH correction is applied'
258      IF ( nn_msshc == 2 ) WRITE(numout,*) '           User defined MSSH correction' 
[2128]259
[3294]260      CALL wrk_dealloc( jpi,jpj, zpromsk )
[2715]261      !
[2128]262   END SUBROUTINE obs_offset_mdt
263 
[2715]264   !!======================================================================
[2128]265END MODULE obs_readmdt
Note: See TracBrowser for help on using the repository browser.