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.
obsinter_z1d.h90 in branches/UKMO/dev_r5518_obs_oper_update_DepthBug/NEMOGCM/NEMO/OPA_SRC/OBS – NEMO

source: branches/UKMO/dev_r5518_obs_oper_update_DepthBug/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90 @ 15365

Last change on this file since 15365 was 15365, checked in by kingr, 12 months ago

Updated logic to catch edge cases.

File size: 8.4 KB
Line 
1   !!----------------------------------------------------------------------
2   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
3   !! $Id$
4   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
5   !!----------------------------------------------------------------------
6
7   SUBROUTINE obs_int_z1d( kpk, kkco, k1dint, kdep, &
8      &                    pobsdep, pobsk, pobs2k,  &
9      &                    pobs, pdep, pobsmask )
10      !!---------------------------------------------------------------------
11      !!
12      !!                   ***  ROUTINE obs_int_z1d ***
13      !!
14      !! ** Purpose : Vertical interpolation to the observation point.
15      !! 
16      !! ** Method  : If k1dint = 0 then use linear interpolation.
17      !!              If k1dint = 1 then use cubic spline interpolation.
18      !!
19      !! ** Action  :
20      !!
21      !! References :
22      !!
23      !! History
24      !!      ! 97-11 (A. Weaver, S. Ricci, N. Daget)
25      !!      ! 06-03 (G. Smith) Conversion to F90 for use with NEMOVAR
26      !!      ! 06-10 (A. Weaver) Cleanup
27      !!      ! 07-01 (K. Mogensen) Use profile rather than single level
28      !!---------------------------------------------------------------------
29
30      !! * Arguments
31      INTEGER, INTENT(IN) :: kpk        ! Number of vertical levels
32      INTEGER, INTENT(IN) :: k1dint     ! 0 = linear; 1 = cubic spline interpolation
33      INTEGER, INTENT(IN) :: kdep       ! Number of levels in profile
34      INTEGER, INTENT(IN), DIMENSION(kdep) :: &
35         & kkco                 ! Array indicies for interpolation
36      REAL(KIND=wp), INTENT(IN), DIMENSION(kdep) :: &
37         & pobsdep              ! Depth of the observation
38      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
39         & pobsk,  &            ! Model profile at a given (lon,lat)
40         & pobs2k, &            ! 2nd derivative of the interpolating function
41         & pdep,   &            ! Model depth array
42         & pobsmask             ! Vertical mask
43      REAL(KIND=wp), INTENT(OUT), DIMENSION(kdep) :: &
44         & pobs                 ! Model equivalent at observation point
45 
46      !! * Local declarations
47      REAL(KIND=wp) :: z1dm       ! Distance above and below obs to model grid points
48      REAL(KIND=wp) :: z1dp         
49      REAL(KIND=wp) :: zsum       ! Dummy variables for computation
50      REAL(KIND=wp) :: zsum2
51      INTEGER :: jdep             ! Observation depths loop variable
52   
53      !------------------------------------------------------------------------
54      ! Loop over all observation depths
55      !------------------------------------------------------------------------
56
57      DO jdep = 1, kdep
58
59         !---------------------------------------------------------------------
60         ! Initialization
61         !---------------------------------------------------------------------
62         z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
63         z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
64         
65         ! Where both levels are masked, return a fill value
66         IF ( pobsmask(kkco(jdep)-1) == 0.0_wp ) .and. (pobsmask(kkco(jdep)) == 0.0_wp) THEN
67            pobs(jdep)  = 99999.
68         ELSE
69         
70            ! Where upper level is masked (e.g., under ice cavity), only use deeper level
71            IF ( pobsmask(kkco(jdep)-1) == 0.0_wp ) z1dm = 0.0_wp
72
73            ! Where deeper level is masked (e.g., near sea bed), only use upper level
74            IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp
75
76            ! Where ob is at or above upper level model T-point, use upper model level
77            ! rather than extrapolate, except where that level is masked
78            IF ( pobsdep(jdep) <= pdep(kkco(jdep)-1) ) .and. &
79            &  ( pobsmask(kkco(jdep)-1) == 0.0_wp ) ) z1dp = 0.0_wp
80
81            ! Where ob is at or below deeper level model T-point, use deeper model level
82            ! rather than extrapolate, except where that level is masked
83            IF ( pobsdep(jdep) >= pdep(kkco(jdep)) ) .and. &
84            &  ( pobsmask(kkco(jdep)) == 0.0_wp ) ) z1dm = 0.0_wp
85
86
87            zsum = z1dm + z1dp
88         
89            IF ( k1dint == 0 ) THEN
90
91               !-----------------------------------------------------------------
92               !  Linear interpolation
93               !-----------------------------------------------------------------
94               pobs(jdep) = (   z1dm * pobsk(kkco(jdep)-1) &
95                  &           + z1dp * pobsk(kkco(jdep)  ) ) / zsum
96
97            ELSEIF ( k1dint == 1 ) THEN
98
99               !-----------------------------------------------------------------
100               ! Cubic spline interpolation
101               !-----------------------------------------------------------------
102               zsum2 = zsum * zsum
103               pobs(jdep)  = (  z1dm                             * pobsk (kkco(jdep)-1) &
104                  &           + z1dp                             * pobsk (kkco(jdep)  ) &
105                  &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco(jdep)-1) &
106                  &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco(jdep)  ) &
107                  &             ) / 6.0_wp                                              &
108                  &          ) / zsum
109
110            ENDIF
111      END DO
112
113   END SUBROUTINE obs_int_z1d
114
115   SUBROUTINE obs_int_z1d_spl( kpk, pobsk, pobs2k, &
116      &                        pdep, pobsmask )
117      !!--------------------------------------------------------------------
118      !!
119      !!                  *** ROUTINE obs_int_z1d_spl ***
120      !!
121      !! ** Purpose : Compute the local vector of vertical second-derivatives
122      !!              of the interpolating function used with a cubic spline.
123      !! 
124      !! ** Method  :
125      !!
126      !!    Top and bottom boundary conditions on the 2nd derivative are
127      !!    set to zero.
128      !!
129      !! ** Action  :
130      !!
131      !! References :
132      !!
133      !! History
134      !!      ! 01-11 (A. Weaver, S. Ricci, N. Daget)
135      !!      ! 06-03 (G. Smith) Conversion to F90 for use with NEMOVAR
136      !!      ! 06-10 (A. Weaver) Cleanup
137      !!----------------------------------------------------------------------
138     
139      !! * Arguments
140      INTEGER, INTENT(IN) :: kpk               ! Number of vertical levels
141      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
142         & pobsk, &          ! Model profile at a given (lon,lat)
143         & pdep,  &          ! Model depth array
144         & pobsmask          ! Vertical mask
145      REAL(KIND=wp), INTENT(OUT), DIMENSION(kpk) :: &
146         & pobs2k            ! 2nd derivative of the interpolating function
147 
148      !! * Local declarations
149      INTEGER :: jk
150      REAL(KIND=wp) :: za
151      REAL(KIND=wp) :: zb
152      REAL(KIND=wp) :: zc
153      REAL(KIND=wp) :: zpa
154      REAL(KIND=wp) :: zkm
155      REAL(KIND=wp) :: zkp
156      REAL(KIND=wp) :: zk
157      REAL(KIND=wp), DIMENSION(kpk-1) :: &
158         & zs, &
159         & zp, &
160         & zu, &
161         & zv
162
163      !-----------------------------------------------------------------------
164      ! Matrix initialisation
165      !-----------------------------------------------------------------------
166      zs(1) =  0.0_wp
167      zp(1) =  0.0_wp
168      zv(1) = -0.5_wp
169      DO jk = 2, kpk-1
170         zs(jk) =  ( pdep(jk  ) - pdep(jk-1) ) &
171            &    / ( pdep(jk+1) - pdep(jk-1) )
172         zp(jk) = zs(jk) * zv(jk-1) + 2.0_wp
173         zv(jk) = ( zs(jk) - 1.0_wp ) / zp(jk)
174      END DO
175 
176      !-----------------------------------------------------------------------
177      ! Solution of the tridiagonal system
178      !-----------------------------------------------------------------------
179 
180      ! Top boundary condition
181      zu(1) = 0.0_wp
182 
183      DO jk = 2, kpk-1
184         za = pdep(jk+1) - pdep(jk-1)
185         zb = pdep(jk+1) - pdep(jk  )
186         zc = pdep(jk  ) - pdep(jk-1)
187 
188         zpa = 6.0_wp / ( zp(jk) * za )
189         zkm = zpa / zc
190         zkp = zpa / zb
191         zk  = - ( zkm + zkp )
192 
193         zu(jk) =  pobsk(jk+1) * zkp  &
194            &    + pobsk(jk  ) * zk   &
195            &    + pobsk(jk-1) * zkm  &
196            &    + zu(jk-1) * ( -zs(jk) / zp(jk) )
197      END DO
198 
199      !-----------------------------------------------------------------------
200      ! Second derivative
201      !-----------------------------------------------------------------------
202      pobs2k(kpk) = 0.0_wp
203 
204      ! Bottom boundary condition
205      DO jk = kpk-1, 1, -1
206         pobs2k(jk) = zv(jk) * pobs2k(jk+1) + zu(jk)
207         IF ( pobsmask(jk+1) == 0.0_wp ) pobs2k(jk) = 0.0_wp
208      END DO
209 
210  END SUBROUTINE obs_int_z1d_spl
211
212
Note: See TracBrowser for help on using the repository browser.