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

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis3/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90 @ 11256

Last change on this file since 11256 was 11256, checked in by rrenshaw, 5 years ago

bug fix for obs below model bathymetry
File size: 8.1 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      !!      ! 19-07 (R. Renshaw) Checks on kkco and zsum to avoid occasional NaNs
29      !!---------------------------------------------------------------------
30
31      !! * Arguments
32      INTEGER, INTENT(IN) :: kpk        ! Number of vertical levels
33      INTEGER, INTENT(IN) :: k1dint     ! 0 = linear; 1 = cubic spline interpolation
34      INTEGER, INTENT(IN) :: kdep       ! Number of levels in profile
35      INTEGER, INTENT(IN), DIMENSION(kdep) :: &
36         & kkco                 ! Array indicies for interpolation
37      REAL(KIND=wp), INTENT(IN), DIMENSION(kdep) :: &
38         & pobsdep              ! Depth of the observation
39      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
40         & pobsk,  &            ! Model profile at a given (lon,lat)
41         & pobs2k, &            ! 2nd derivative of the interpolating function
42         & pdep,   &            ! Model depth array
43         & pobsmask             ! Vertical mask
44      REAL(KIND=wp), INTENT(OUT), DIMENSION(kdep) :: &
45         & pobs                 ! Model equivalent at observation point
46 
47      !! * Local declarations
48      REAL(KIND=wp) :: z1dm       ! Distance above and below obs to model grid points
49      REAL(KIND=wp) :: z1dp         
50      REAL(KIND=wp) :: zsum       ! Dummy variables for computation
51      REAL(KIND=wp) :: zsum2
52      INTEGER :: jdep             ! Observation depths loop variable
53      INTEGER :: kkco_fix        ! kkco(jdep) limited to number of model levels
54   
55      !------------------------------------------------------------------------
56      ! Loop over all observation depths
57      !------------------------------------------------------------------------
58
59      DO jdep = 1, kdep
60
61         !---------------------------------------------------------------------
62         ! Initialization
63         !---------------------------------------------------------------------
64!        kkco (from obs_level_search) sometimes has model level kpk+1 which doesn't exist
65         kkco_fix = MIN( kkco(jdep), kpk )
66         z1dm = ( pdep(kkco_fix) - pobsdep(jdep)      )
67         z1dp = ( pobsdep(jdep)    - pdep(kkco_fix-1) )
68!        Where ob is below model bottom, use model bottom and don't extrapolate
69         IF ( pdep(kkco_fix) <= pobsdep(jdep) ) z1dm = 0.0_wp
70!        Where lower level is missing, use higher level
71         IF ( pobsmask(kkco_fix) == 0.0_wp ) z1dp = 0.0_wp
72
73         zsum = z1dm + z1dp
74         ! if pobsmask==0 and model level depth==observed depth, we get zsum=0
75         IF ( zsum > 0.0_wp ) THEN
76         
77         IF ( k1dint == 0 ) THEN
78
79            !-----------------------------------------------------------------
80            !  Linear interpolation
81            !-----------------------------------------------------------------
82            pobs(jdep) = (   z1dm * pobsk(kkco_fix-1) &
83               &           + z1dp * pobsk(kkco_fix  ) ) / zsum
84
85         ELSEIF ( k1dint == 1 ) THEN
86
87            !-----------------------------------------------------------------
88            ! Cubic spline interpolation
89            !-----------------------------------------------------------------
90            zsum2 = zsum * zsum
91            pobs(jdep)  = (  z1dm                             * pobsk (kkco_fix-1) &
92               &           + z1dp                             * pobsk (kkco_fix  ) &
93               &           + ( z1dm * ( z1dm * z1dm - zsum2 ) * pobs2k(kkco_fix-1) &
94               &           +   z1dp * ( z1dp * z1dp - zsum2 ) * pobs2k(kkco_fix  ) &
95               &             ) / 6.0_wp                                              &
96               &          ) / zsum
97
98         ENDIF
99
100         ELSE  ! take value directly from the higher model level
101            pobs(jdep)  = pobsk(kkco_fix-1)
102         ENDIF
103
104      END DO
105
106   END SUBROUTINE obs_int_z1d
107
108   SUBROUTINE obs_int_z1d_spl( kpk, pobsk, pobs2k, &
109      &                        pdep, pobsmask )
110      !!--------------------------------------------------------------------
111      !!
112      !!                  *** ROUTINE obs_int_z1d_spl ***
113      !!
114      !! ** Purpose : Compute the local vector of vertical second-derivatives
115      !!              of the interpolating function used with a cubic spline.
116      !! 
117      !! ** Method  :
118      !!
119      !!    Top and bottom boundary conditions on the 2nd derivative are
120      !!    set to zero.
121      !!
122      !! ** Action  :
123      !!
124      !! References :
125      !!
126      !! History
127      !!      ! 01-11 (A. Weaver, S. Ricci, N. Daget)
128      !!      ! 06-03 (G. Smith) Conversion to F90 for use with NEMOVAR
129      !!      ! 06-10 (A. Weaver) Cleanup
130      !!----------------------------------------------------------------------
131     
132      !! * Arguments
133      INTEGER, INTENT(IN) :: kpk               ! Number of vertical levels
134      REAL(KIND=wp), INTENT(IN), DIMENSION(kpk) :: &
135         & pobsk, &          ! Model profile at a given (lon,lat)
136         & pdep,  &          ! Model depth array
137         & pobsmask          ! Vertical mask
138      REAL(KIND=wp), INTENT(OUT), DIMENSION(kpk) :: &
139         & pobs2k            ! 2nd derivative of the interpolating function
140 
141      !! * Local declarations
142      INTEGER :: jk
143      REAL(KIND=wp) :: za
144      REAL(KIND=wp) :: zb
145      REAL(KIND=wp) :: zc
146      REAL(KIND=wp) :: zpa
147      REAL(KIND=wp) :: zkm
148      REAL(KIND=wp) :: zkp
149      REAL(KIND=wp) :: zk
150      REAL(KIND=wp), DIMENSION(kpk-1) :: &
151         & zs, &
152         & zp, &
153         & zu, &
154         & zv
155
156      !-----------------------------------------------------------------------
157      ! Matrix initialisation
158      !-----------------------------------------------------------------------
159      zs(1) =  0.0_wp
160      zp(1) =  0.0_wp
161      zv(1) = -0.5_wp
162      DO jk = 2, kpk-1
163         zs(jk) =  ( pdep(jk  ) - pdep(jk-1) ) &
164            &    / ( pdep(jk+1) - pdep(jk-1) )
165         zp(jk) = zs(jk) * zv(jk-1) + 2.0_wp
166         zv(jk) = ( zs(jk) - 1.0_wp ) / zp(jk)
167      END DO
168 
169      !-----------------------------------------------------------------------
170      ! Solution of the tridiagonal system
171      !-----------------------------------------------------------------------
172 
173      ! Top boundary condition
174      zu(1) = 0.0_wp
175 
176      DO jk = 2, kpk-1
177         za = pdep(jk+1) - pdep(jk-1)
178         zb = pdep(jk+1) - pdep(jk  )
179         zc = pdep(jk  ) - pdep(jk-1)
180 
181         zpa = 6.0_wp / ( zp(jk) * za )
182         zkm = zpa / zc
183         zkp = zpa / zb
184         zk  = - ( zkm + zkp )
185 
186         zu(jk) =  pobsk(jk+1) * zkp  &
187            &    + pobsk(jk  ) * zk   &
188            &    + pobsk(jk-1) * zkm  &
189            &    + zu(jk-1) * ( -zs(jk) / zp(jk) )
190      END DO
191 
192      !-----------------------------------------------------------------------
193      ! Second derivative
194      !-----------------------------------------------------------------------
195      pobs2k(kpk) = 0.0_wp
196 
197      ! Bottom boundary condition
198      DO jk = kpk-1, 1, -1
199         pobs2k(jk) = zv(jk) * pobs2k(jk+1) + zu(jk)
200         IF ( pobsmask(jk+1) == 0.0_wp ) pobs2k(jk) = 0.0_wp
201      END DO
202 
203  END SUBROUTINE obs_int_z1d_spl
204
205
Note: See TracBrowser for help on using the repository browser.