1 | subroutine interpolate(ic, jc, kc, im, jm, km, fx,fy,fz, E_k,nft) |
---|
2 | |
---|
3 | use module3d_phy |
---|
4 | use tracer_vars |
---|
5 | !cdc use climat_perturb_mois_mod ! on en a besoin, notamment pour NFT |
---|
6 | |
---|
7 | implicit none |
---|
8 | integer, intent(IN) :: ic, jc, kc, im, jm, km |
---|
9 | real, intent(INOUT) :: fx, fy, fz, E_k |
---|
10 | real :: fx_1, fy_1, fz_1, a_xy, a_xy_p, f_zt, a_tmp, ac , al |
---|
11 | real :: fra1, fra2, fra0, accucumul1, accucumul2, accucumul0, & |
---|
12 | thinning, thinning2, f_k, a_xy_m, fzz, fprime, fsec |
---|
13 | integer :: index1, index2, index0 |
---|
14 | integer :: nft |
---|
15 | |
---|
16 | integer :: kk ! indice vertical pour definition de E, mais on veut conserver la valeur de k |
---|
17 | ! real,dimension(nz) :: E ! vertical coordinate in ice, scaled to H zeta |
---|
18 | |
---|
19 | !! Note that input (im,jm,km) refer to points on (ij,jj) grid |
---|
20 | !! whereas CellTest search was performed on (i,j) grid |
---|
21 | |
---|
22 | |
---|
23 | E(1)=0. |
---|
24 | E(NZ)=1. |
---|
25 | do KK=1,NZ |
---|
26 | if ((KK.ne.1).and.(KK.ne.NZ)) E(KK)=(KK-1.)/(NZ-1.) |
---|
27 | end do |
---|
28 | |
---|
29 | fx_1 = 1.0-fx |
---|
30 | fy_1 = 1.0-fy |
---|
31 | fz_1 = 1.0-fz |
---|
32 | |
---|
33 | ! start with bilinear interp in horizontal plane |
---|
34 | ! above |
---|
35 | a_xy = fx_1*fy_1*tdepk(km,im,jm)+fx*fy_1*tdepk(km,im+1,jm)+ & |
---|
36 | fx_1*fy*tdepk(km,im,jm+1)+fx*fy*tdepk(km,im+1,jm+1) |
---|
37 | ! below |
---|
38 | a_xy_p = fx_1*fy_1*tdepk(km+1,im,jm)+fx*fy_1*tdepk(km+1,im+1,jm)+ & |
---|
39 | fx_1*fy*tdepk(km+1,im,jm+1)+fx*fy*tdepk(km+1,im+1,jm+1) |
---|
40 | |
---|
41 | |
---|
42 | if ( ((a_xy - a_xy_p) < 200.).or.(H(ic+1,jc+1)<100.)) then |
---|
43 | tdep(kc,ic,jc) = fz_1*a_xy+fz*a_xy_p |
---|
44 | |
---|
45 | else if ((a_xy_p<time_max_accu +1000.) & |
---|
46 | .or.(fz<1.E-7).or.(fz_1<1.E-7) ) then ! added 25/03/04: don t bother for < 1/50mm |
---|
47 | |
---|
48 | a_tmp = fz_1*a_xy+fz*a_xy_p |
---|
49 | |
---|
50 | if ((km==1).and.(fz<0.5)) then |
---|
51 | tdep(kc,ic,jc) = a_tmp |
---|
52 | |
---|
53 | else ! cubic interp from the closest neighbour |
---|
54 | if (fz>0.5) then ! furthest to km, use km+2 unless bottom |
---|
55 | a_xy_m= a_xy |
---|
56 | a_xy = a_xy_p |
---|
57 | if (km+1==nz) then |
---|
58 | a_xy_p = a_xy + 1.3*(a_xy - a_xy_m) ! like articially old 22nd layer |
---|
59 | else |
---|
60 | a_xy_p = fx_1*fy_1*tdepk(km+2,im,jm)+fx*fy_1*tdepk(km+2,im+1,jm)+ & |
---|
61 | fx_1*fy*tdepk(km+2,im,jm+1)+fx*fy*tdepk(km+2,im+1,jm+1) |
---|
62 | end if |
---|
63 | fzz= - fz_1 ! because km+1 closer than km and E_D<E(km+1) |
---|
64 | else |
---|
65 | a_xy_m = fx_1*fy_1*tdepk(km-1,im,jm)+fx*fy_1*tdepk(km-1,im+1,jm)+ & |
---|
66 | fx_1*fy*tdepk(km-1,im,jm+1)+fx*fy*tdepk(km-1,im+1,jm+1) |
---|
67 | fzz = fz |
---|
68 | end if |
---|
69 | fprime = - abs(a_xy_p - a_xy_m)/2. ! /E(km+1)-E(km) included in fz |
---|
70 | fsec = - abs(a_xy_p - 2.*a_xy + a_xy_m) ! idem |
---|
71 | ! there is a risk of overshooting, limit change to 1.3x linear change |
---|
72 | al = abs(tdepk(kc,ic,jc) - a_tmp) |
---|
73 | ac = abs(tdepk(kc,ic,jc) - (a_xy + fzz*fprime + fzz**2*fsec/2.)) |
---|
74 | if (ac<(1.3*al)) then |
---|
75 | tdep(kc,ic,jc) = a_xy + fzz*fprime + fzz**2*fsec/2. |
---|
76 | else |
---|
77 | tdep(kc,ic,jc) = a_tmp |
---|
78 | end if |
---|
79 | end if |
---|
80 | |
---|
81 | |
---|
82 | else |
---|
83 | |
---|
84 | index1 = floor(-a_xy/100.) ! gives the youngest: ind1+1 < tdepk < ind1 |
---|
85 | fra1 = real((-a_xy/100.) - index1) |
---|
86 | accucumul1 = real( -fra1*(accucumul(index1)-accucumul(index1+1)) + accucumul(index1) ) |
---|
87 | index2 = floor(-a_xy_p/100.) |
---|
88 | fra2 = real((-a_xy_p/100.) - index2) |
---|
89 | accucumul2 = real( -fra2 * (accucumul(index2)-accucumul(index2+1)) + accucumul(index2) ) |
---|
90 | |
---|
91 | !print*, 'ic,jc,kc :',ic,jc,kc |
---|
92 | !print*, 'a_xy et a_xy_p', a_xy, a_xy_p |
---|
93 | !print*, 'index1 et 2',index1,index2 |
---|
94 | !print*, 'fra1 et 2', fra1,fra2 |
---|
95 | !print*, 'accucumul 1 et 2 :', accucumul1, accucumul2,accucumul2-accucumul1 !aurel neem |
---|
96 | thinning = ( real( -(E(km+1)-E(km)) / (accucumul2 - accucumul1)) ) ! checked |
---|
97 | !print*, 'thinning ', thinning !aurel neem |
---|
98 | |
---|
99 | if (fz < 0.5) then |
---|
100 | !======= closer to km, probably upward motion, rather seldom NL thinks |
---|
101 | if (km>1) then |
---|
102 | a_xy_m = fx_1*fy_1*tdepk(km-1,im,jm)+fx*fy_1*tdepk(km-1,im+1,jm)+ & |
---|
103 | fx_1*fy*tdepk(km-1,im,jm+1)+fx*fy*tdepk(km-1,im+1,jm+1) |
---|
104 | ! should be a_xy_m - a_xy > 0 |
---|
105 | if ((a_xy_m>time_max_accu+ 500.).and.( (a_xy_m - a_xy) > 200.)) then |
---|
106 | index0 = floor(-a_xy_m/100.) |
---|
107 | fra0 = real((-a_xy_m/100.) - index0) |
---|
108 | accucumul0 = real( -fra0 * (accucumul(index0)-accucumul(index0+1)) & |
---|
109 | + accucumul(index0) ) |
---|
110 | ! acc0 > acc1 |
---|
111 | |
---|
112 | !print*, 'ic,jc,kc :',ic,jc,kc |
---|
113 | !print*, 'accucumul 0 et 1 : ', accucumul0, accucumul1,accucumul1-accucumul0 !aurel neem |
---|
114 | thinning2 = ( real( -(E(km)-E(km-1)) / (accucumul1 - accucumul0)) ) ! ok |
---|
115 | !print*, 'thinning2 ', thinning2 !aurel neem |
---|
116 | |
---|
117 | thinning = max(1.e-8, (0.5-fz/2.0)*thinning2 + (0.5+fz/2.0)*thinning) |
---|
118 | |
---|
119 | ! fz/c2 says my theory, though for Greenland I find that /c3 or /c4 may be better... |
---|
120 | end if |
---|
121 | end if |
---|
122 | |
---|
123 | f_k = real( accucumul1 - (E_k - E(km) ) / thinning ) ! ok: f_k < acc1 |
---|
124 | ! downward search |
---|
125 | index0=index1-1 ! no inversion in stratigraphy |
---|
126 | L4: do |
---|
127 | index0 = index0 + 1 |
---|
128 | if ((index0>=NFT-1).or.(accucumul(index0) < f_k )) exit L4 |
---|
129 | ! index < ind(t) < index-1, t(i+1) < t < t(i) |
---|
130 | end do L4 |
---|
131 | index0 = index0-1 |
---|
132 | |
---|
133 | !========================================= |
---|
134 | else ! fz> 0.5, probably downwards |
---|
135 | |
---|
136 | if (km<nz-1) then |
---|
137 | a_xy_m = fx_1*fy_1*tdepk(km+2,im,jm)+fx*fy_1*tdepk(km+2,im+1,jm)+ & |
---|
138 | fx_1*fy*tdepk(km+2,im,jm+1)+fx*fy*tdepk(km+2,im+1,jm+1) |
---|
139 | if ((a_xy_m>time_max_accu+ 500.).and.( (a_xy_m - a_xy_p) < -200.)) then |
---|
140 | index0 = floor(-a_xy_m/100.) |
---|
141 | fra0 = real((-a_xy_m/100.) - index0) |
---|
142 | accucumul0 = real( -fra0 * (accucumul(index0)-accucumul(index0+1)) & |
---|
143 | + accucumul(index0) ) |
---|
144 | ! print*, 'ic,jc,kc :',ic,jc,kc |
---|
145 | ! print*, 'accucumul 0 et 2 : ', accucumul0, accucumul2,accucumul0-accucumul2 !aurel neem |
---|
146 | thinning2 = ( real( -(E(km+2)-E(km+1)) / (accucumul0 - accucumul2)) ) |
---|
147 | ! print*, 'thinning2 ', thinning2 !aurel neem |
---|
148 | ! E decreases upwards! |
---|
149 | thinning = max(1.e-8, (0.5-fz_1/2.0)*thinning2 + (0.5+fz_1/2.0)*thinning) |
---|
150 | end if |
---|
151 | end if |
---|
152 | |
---|
153 | f_k = real( accucumul2 + (E(km+1) - E_k) / thinning ) |
---|
154 | ! upward search |
---|
155 | index0=index2+2 ! no inversion in stratigraphy |
---|
156 | L5: do |
---|
157 | index0 = index0 - 1 |
---|
158 | if ((index0<=1).or.(accucumul(index0)> f_k )) exit L5 |
---|
159 | ! index < ind(t) < index+1, t(i+1) < t < t(i) |
---|
160 | end do L5 |
---|
161 | end if |
---|
162 | ! now we know that we have index-1 < tdep(k) < index from accucumul condition |
---|
163 | |
---|
164 | if ((index0==1).or.(index0==NFT-1)) then |
---|
165 | fzz=2. |
---|
166 | else |
---|
167 | fzz = real( (accucumul(index0) - f_k)/ (accucumul(index0)-accucumul(index0+1)) ) |
---|
168 | end if |
---|
169 | |
---|
170 | if ((fzz>1.0).or.(fzz<0.0)) then |
---|
171 | tdep(kc,ic,jc) = fz_1*a_xy + fz*a_xy_p |
---|
172 | else |
---|
173 | tdep(kc,ic,jc) = ( fzz + real(index0) ) * real(-100) |
---|
174 | end if |
---|
175 | |
---|
176 | end if |
---|
177 | |
---|
178 | !=============== end tdep-interp ==================== |
---|
179 | |
---|
180 | ! f_zt = (tdep(kc,ic,jc) - a_xy) / (a_xy_p - a_xy) |
---|
181 | ! NL 21/11/02: use only fz, not f_zt, i don't think f_zt is right |
---|
182 | |
---|
183 | a_xy = fx_1*fy_1*xdepk(km,im,jm)+fx*fy_1*xdepk(km,im+1,jm)+ & |
---|
184 | fx_1*fy*xdepk(km,im,jm+1)+fx*fy*xdepk(km,im+1,jm+1) |
---|
185 | a_xy_p = fx_1*fy_1*xdepk(km+1,im,jm)+fx*fy_1*xdepk(km+1,im+1,jm)+ & |
---|
186 | fx_1*fy*xdepk(km+1,im,jm+1)+fx*fy*xdepk(km+1,im+1,jm+1) |
---|
187 | |
---|
188 | xdep(kc,ic,jc) = fz_1*a_xy+fz*a_xy_p |
---|
189 | |
---|
190 | !========= YDEP =================== |
---|
191 | |
---|
192 | a_xy = fx_1*fy_1*ydepk(km,im,jm)+fx*fy_1*ydepk(km,im+1,jm)+ & |
---|
193 | fx_1*fy*ydepk(km,im,jm+1)+fx*fy*ydepk(km,im+1,jm+1) |
---|
194 | a_xy_p = fx_1*fy_1*ydepk(km+1,im,jm)+fx*fy_1*ydepk(km+1,im+1,jm)+ & |
---|
195 | fx_1*fy*ydepk(km+1,im,jm+1)+fx*fy*ydepk(km+1,im+1,jm+1) |
---|
196 | |
---|
197 | ydep(kc,ic,jc) = fz_1*a_xy+fz*a_xy_p |
---|
198 | |
---|
199 | |
---|
200 | end subroutine interpolate |
---|
201 | !========================================================================== |
---|