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