source: branches/GRISLIv3/SOURCES/interpolate_tracer.f90 @ 409

Last change on this file since 409 was 409, checked in by dumas, 15 months ago

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