source: utils/tools/ABL_TOOLS/module_interp.F90 @ 11589

Last change on this file since 11589 was 11589, checked in by gsamson, 12 months ago

dev_r11265_ABL : see #2131

  • ABL_TOOLS first working version (README empty and arch files ignored for now)
File size: 18.3 KB
Line 
1MODULE module_interp
2   !!======================================================================
3   !!                   ***  MODULE  module_interp  ***
4   !! Ocean forcing:  bulk thermohaline forcing of the ocean (or ice)
5   !!=====================================================================
6   !! History : 2016-10  (F. Lemarié)  Original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   zinterp                :
11   !!   reconstructandremap    :
12   !!   reconstructandremap_ps :
13   !!
14   !!----------------------------------------------------------------------
15   IMPLICIT NONE
16
17CONTAINS
18
19   SUBROUTINE zinterp( jpi, jpj, jpka, jpka_in, ind, tab_in, e3t_in, e3_bak, &   
20      &                              e3t_out, tab_out, interp_type           )   
21     
22      !!---------------------------------------------------------------------
23      !!                    ***  ROUTINE zinterp  ***
24      !!                   
25      !! ** Purpose :   
26      !!
27      !! ** Method  :
28      !!
29      !! ** Action  : 
30      !!----------------------------------------------------------------------
31      INTEGER, INTENT(in   ) ::   jpi, jpj
32      INTEGER, INTENT(in   ) ::   jpka, jpka_in, interp_type
33      INTEGER, INTENT(in   ) ::   ind     ( 1:jpi, 1:jpj            )
34      REAL(8), INTENT(inout) ::   tab_in  ( 1:jpi, 1:jpj, 1:jpka_in )
35      REAL(8), INTENT(in   ) ::   e3t_in  ( 1:jpi, 1:jpj, 1:jpka_in )
36      REAL(8), INTENT(in   ) ::   e3_bak  ( 1:jpi, 1:jpj            )           
37      REAL(8), INTENT(in   ) ::   e3t_out (               1:jpka+1  )                           
38      REAL(8), INTENT(  out) ::   tab_out ( 1:jpi, 1:jpj, 1:jpka+1  ) 
39      !!
40      INTEGER                :: ji,jj,k_in
41      REAL(8)                :: val1,val2,cff
42     
43      SELECT CASE(interp_type)
44      CASE(1)      ! WENO
45         DO jj = 1,jpj
46            DO ji = 1,jpi
47               k_in = ind( ji, jj )
48               val1 = tab_in ( ji, jj, k_in-1 )
49               val2 = tab_in ( ji, jj, k_in   )
50               cff  =        val1 * e3_bak( ji, jj         )     &
51                  & +        val2 * e3t_in( ji, jj, k_in-1 )     &
52                  & + (val2-val1) * e3t_in( ji, jj, k_in   )             
53               tab_in( ji, jj, k_in ) = cff / ( e3_bak( ji, jj ) + e3t_in ( ji, jj, k_in-1 ) )
54               !
55               CALL reconstructandremap( tab_in ( ji, jj, 1:k_in   ), e3t_in( ji, jj, 1:k_in ),   &
56                  &                      tab_out( ji, jj, 2:jpka+1 ), e3t_out (    2:jpka+1  ),   &
57                  &                      k_in, jpka ) 
58               !
59            END DO
60         END DO   
61      CASE(2)      ! SPLINES
62         DO jj = 1,jpj
63            DO ji = 1,jpi
64               k_in = ind( ji, jj )
65               val1 = tab_in ( ji, jj, k_in-1 )
66               val2 = tab_in ( ji, jj, k_in   )
67               cff  =        val1 * e3_bak( ji, jj         )     &
68                  & +        val2 * e3t_in( ji, jj, k_in-1 )     &
69                  & + (val2-val1) * e3t_in( ji, jj, k_in   )             
70               tab_in( ji, jj, k_in ) = cff / ( e3_bak( ji, jj ) + e3t_in ( ji, jj, k_in-1 ) )
71               !
72               CALL reconstructandremap_ps( tab_in ( ji, jj, 1:k_in   ), e3t_in( ji, jj, 1:k_in ),   &
73                  &                         tab_out( ji, jj, 2:jpka+1 ), e3t_out (    2:jpka+1  ),   &
74                  &                         k_in, jpka ) 
75               !
76            END DO
77         END DO 
78      CASE DEFAULT
79         WRITE(*,*) "### Error: problem in zinterp, interp_type not set properly"
80         STOP
81      END SELECT     
82      !
83   END SUBROUTINE zinterp     
84
85
86
87
88
89!
90!===================================================================================================
91subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout)
92!---------------------------------------------------------------------------------------------------
93      implicit none
94      integer            :: N, Nout
95      real(8)            :: tabin(N), tabout(Nout)
96      real(8)            :: hin(N), hout(Nout)
97      real(8)            :: coeffremap(N,3),zwork(N,3)
98      real(8)            :: zwork2(N+1,3)
99      integer            :: k
100      real(8), parameter :: dsmll=1.0d-8 
101      real(8)            :: q,q01,q02,q001,q002,q0
102      real(8)            :: z_win(1:N+1), z_wout(1:Nout+1)
103      real(8),parameter  :: dpthin = 1.D-3
104      integer            :: k1, kbox, ktop, ka, kbot
105      real(8)            :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
106!-----
107
108!----------------------
109      z_win(1)=0.; z_wout(1)= 0.
110      do k=1,N
111       z_win(k+1)=z_win(k)+hin(k)
112      enddo 
113     
114      do k=1,Nout
115       z_wout(k+1)=z_wout(k)+hout(k)       
116      enddo       
117
118        do k=2,N
119          zwork(k,1)=1./(hin(k-1)+hin(k))
120        enddo
121       
122        do k=2,N-1
123          q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
124          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
125          zwork(k,3)=q0
126        enddo       
127     
128        do k= 2,N
129        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
130        enddo
131       
132        coeffremap(:,1) = tabin(:)
133 
134         do k=2,N-1
135        q001 = hin(k)*zwork2(k+1,1)
136        q002 = hin(k)*zwork2(k,1)       
137        if (q001*q002 < 0) then
138          q001 = 0.
139          q002 = 0.
140        endif
141        q=zwork(k,2)
142        q01=q*zwork2(k+1,1)
143        q02=q*zwork2(k,1)
144        if (abs(q001) > abs(q02)) q001 = q02
145        if (abs(q002) > abs(q01)) q002 = q01
146       
147        q=(q001-q002)*zwork(k,3)
148        q001=q001-q*hin(k+1)
149        q002=q002+q*hin(k-1)
150       
151        coeffremap(k,3)=coeffremap(k,1)+q001
152        coeffremap(k,2)=coeffremap(k,1)-q002
153       
154        zwork2(k,1)=(2.*q001-q002)**2
155        zwork2(k,2)=(2.*q002-q001)**2
156        enddo
157       
158        do k=1,N
159        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then
160        coeffremap(k,3) = coeffremap(k,1)
161        coeffremap(k,2) = coeffremap(k,1)
162        zwork2(k,1) = 0.
163        zwork2(k,2) = 0.
164        endif
165        enddo
166       
167        do k=2,N
168        q002=max(zwork2(k-1,2),dsmll)
169        q001=max(zwork2(k,1),dsmll)
170        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002)
171        enddo
172       
173        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
174        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
175 
176        do k=1,N
177        q01=zwork2(k+1,3)-coeffremap(k,1)
178        q02=coeffremap(k,1)-zwork2(k,3)
179        q001=2.*q01
180        q002=2.*q02
181        if (q01*q02<0) then
182          q01=0.
183          q02=0.
184        elseif (abs(q01)>abs(q002)) then
185          q01=q002
186        elseif (abs(q02)>abs(q001)) then
187          q02=q001
188        endif
189        coeffremap(k,2)=coeffremap(k,1)-q02
190        coeffremap(k,3)=coeffremap(k,1)+q01
191        enddo
192
193      zbot=0.0
194      kbot=1
195      do k=1,Nout
196        ztop=zbot  !top is bottom of previous layer
197        ktop=kbot
198        if     (ztop.ge.z_win(ktop+1)) then
199          ktop=ktop+1
200        endif
201       
202        zbot=z_wout(k+1)
203        zthk=zbot-ztop
204
205        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then
206
207          kbot=ktop
208          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
209            kbot=kbot+1
210          enddo
211          zbox=zbot
212          do k1= k+1,Nout
213            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then
214              exit !thick layer
215            else
216              zbox=z_wout(k1+1)  !include thin adjacent layers
217              if     (zbox.eq.z_wout(Nout+1)) then
218                exit !at bottom
219              endif
220            endif
221          enddo
222          zthk=zbox-ztop
223
224          kbox=ktop
225          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
226            kbox=kbox+1
227          enddo
228         
229          if     (ktop.eq.kbox) then
230
231
232            if     (z_wout(k)  .ne.z_win(kbox)   .or.z_wout(k+1).ne.z_win(kbox+1)     ) then
233
234              if     (hin(kbox).gt.dpthin) then
235                q001 = (zbox-z_win(kbox))/hin(kbox)
236                q002 = (ztop-z_win(kbox))/hin(kbox)
237                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
238                q02=q01-1.+(q001+q002)
239                q0=1.-q01-q02
240              else
241                q0 = 1.0
242                q01 = 0.
243                q02 = 0.
244              endif
245          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
246             
247            else
248            tabout(k) = tabin(kbox)
249             
250            endif
251
252          else
253
254            if     (ktop.le.k .and. kbox.ge.k) then
255              ka = k
256            elseif (kbox-ktop.ge.3) then
257              ka = (kbox+ktop)/2
258            elseif (hin(ktop).ge.hin(kbox)) then
259              ka = ktop
260            else
261              ka = kbox
262            endif !choose ka
263
264            offset=coeffremap(ka,1)
265
266            qtop = z_win(ktop+1)-ztop !partial layer thickness
267            if     (hin(ktop).gt.dpthin) then
268              q=(ztop-z_win(ktop))/hin(ktop)
269              q01=q*(q-1.)
270              q02=q01+q
271              q0=1-q01-q02           
272            else
273              q0 = 1.
274              q01 = 0.
275              q02 = 0.
276            endif
277           
278            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
279
280            do k1= ktop+1,kbox-1
281              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
282            enddo !k1
283
284           
285            qbot = zbox-z_win(kbox) !partial layer thickness
286            if     (hin(kbox).gt.dpthin) then
287              q=qbot/hin(kbox)
288              q01=(q-1.)**2
289              q02=q01-1.+q
290              q0=1-q01-q02                           
291            else
292              q0 = 1.0
293              q01 = 0.
294              q02 = 0.
295            endif
296           
297            tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
298           
299            rpsum=1.0d0/zthk
300              tabout(k)=offset+tsum*rpsum
301             
302          endif !single or multiple layers
303        else
304        if (k==1) then
305        print *,'problem = ',zthk,z_wout(k+1),hout(1)
306        endif
307         tabout(k) = tabout(k-1)
308         
309        endif !normal:thin layer
310      enddo !k
311           
312      return
313     
314!---------------------------------------------------------------------------------------------------
315end subroutine reconstructandremap
316!===================================================================================================
317!
318
319
320
321
322
323
324
325
326
327!
328!===================================================================================================
329      subroutine reconstructandremap_ps(tabin,hin,tabout,hout,N,Nout)      ! parabloc spline
330!---------------------------------------------------------------------------------------------------     
331      implicit none
332      integer N, Nout
333      real(8) tabin(N), tabout(Nout)
334      real(8) hin(N), hout(Nout)
335      real(8) coeffremap(N,3),zwork(N,3)
336      real(8) zwork2(N+1,3)
337     
338      real(8) my_zwork(0:N,3)
339      real(8) my_zwork2(0:N,3)
340     
341      integer k
342      double precision, parameter :: dsmll=1.0d-8 
343      real(8) q,q01,q02,q001,q002,q0
344      real(8) z_win(1:N+1), z_wout(1:Nout+1)
345      real(8),parameter :: dpthin = 1.D-3
346      integer :: k1, kbox, ktop, ka, kbot
347      real(8) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
348      real(8) :: p
349      real(8) :: qtri(0:N)
350
351      z_win(1)=0.; z_wout(1)= 0.
352      do k=1,N
353       z_win(k+1)=z_win(k)+hin(k)
354      enddo 
355     
356      do k=1,Nout
357       z_wout(k+1)=z_wout(k)+hout(k)       
358      enddo       
359
360        do k=2,N
361          zwork(k,1)=1./(hin(k-1)+hin(k))
362        enddo
363       
364        do k=2,N-1
365          q0 = 1./(hin(k-1)+hin(k)+hin(k+1))
366          zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1)
367          zwork(k,3)=q0
368        enddo       
369     
370        do k= 2,N
371        zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1))
372        enddo
373       
374        coeffremap(:,1) = tabin(:)
375 
376         do k=2,N-1
377        q001 = hin(k)*zwork2(k+1,1)
378        q002 = hin(k)*zwork2(k,1)       
379 !       if (q001*q002 < 0) then
380 !         q001 = 0.
381 !         q002 = 0.
382 !       endif
383        q=zwork(k,2)
384        q01=q*zwork2(k+1,1)
385        q02=q*zwork2(k,1)
386 !       if (abs(q001) > abs(q02)) q001 = q02
387 !       if (abs(q002) > abs(q01)) q002 = q01
388       
389        q=(q001-q002)*zwork(k,3)
390        q001=q001-q*hin(k+1)
391        q002=q002+q*hin(k-1)
392       
393        coeffremap(k,3)=coeffremap(k,1)+q001
394        coeffremap(k,2)=coeffremap(k,1)-q002
395       
396        zwork2(k,1)=(2.*q001-q002)**2
397        zwork2(k,2)=(2.*q002-q001)**2
398        enddo
399       
400        do k=1,N
401        if     (k.eq.1 .or. k.eq.N .or.   hin(k).le.dpthin) then
402        coeffremap(k,3) = coeffremap(k,1)
403        coeffremap(k,2) = coeffremap(k,1)
404        zwork2(k,1) = 0.
405        zwork2(k,2) = 0.
406        endif
407        enddo
408       
409        do k=2,N
410        q002=max(zwork2(k-1,2),dsmll)
411        q001=max(zwork2(k,1),dsmll)
412        zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002)
413        enddo
414       
415        zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
416        zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
417 
418        do k=1,N
419        q01=zwork2(k+1,3)-coeffremap(k,1)
420        q02=coeffremap(k,1)-zwork2(k,3)
421!        q001=2.*q01
422!        q002=2.*q02
423!        if (q01*q02<0) then
424!          q01=0.
425!          q02=0.
426!        elseif (abs(q01)>abs(q002)) then
427!          q01=q002
428!        elseif (abs(q02)>abs(q001)) then
429!          q02=q001
430!        endif
431        coeffremap(k,2)=coeffremap(k,1)-q02
432        coeffremap(k,3)=coeffremap(k,1)+q01
433        enddo
434       
435       
436        do k=0,N
437        if (k==0) then
438        my_zwork(k,1)=0.
439        my_zwork(k,2)=1.
440        my_zwork(k,3)=0.5
441        my_zwork2(k,1)=1.5*tabin(1)
442        elseif (k==N) then
443        my_zwork(k,1)=0.5
444        my_zwork(k,2)=1.
445        my_zwork(k,3)=0.
446        my_zwork2(k,1)=1.5*tabin(k)       
447        else
448        my_zwork(k,1)=hin(k+1)
449        my_zwork(k,2)=2.*(hin(k)+hin(k+1))
450        my_zwork(k,3)=hin(k)
451        my_zwork2(k,1)=3.*(hin(k+1)*tabin(k)+hin(k)*tabin(k+1))
452        my_zwork2(k,2)=my_zwork2(k,1)
453        endif
454        enddo
455       
456        qtri(0)=-my_zwork(0,3)/my_zwork(0,2)
457        my_zwork2(0,1)=my_zwork2(0,1)/my_zwork(0,2)
458       
459        do k=1,N
460        p=1.0/(my_zwork(k,2)+my_zwork(k,1)*qtri(k-1))
461        qtri(k)=-my_zwork(k,3)*p
462        my_zwork2(k,1)=(my_zwork2(k,1)-my_zwork(k,1)*my_zwork2(k-1,1))*p
463        enddo
464       
465        do k=N-1,0,-1
466        my_zwork2(k,1)=my_zwork2(k,1)+qtri(k)*my_zwork2(k+1,1)
467        enddo
468       
469        do k=1,N
470        coeffremap(k,2)=my_zwork2(k-1,1)
471        coeffremap(k,3)=my_zwork2(k,1)
472        enddo
473       
474        do k=2,N-1
475!        print *,'VAL22 = ',my_zwork(k,1)*my_zwork2(k-1,1)
476!     &+my_zwork(k,2)*my_zwork2(k,1)
477!     &  +my_zwork(k,3)*my_zwork2(k+1,1),my_zwork2(k,2)
478        enddo
479
480      zbot=0.0
481      kbot=1
482      do k=1,Nout
483        ztop=zbot  !top is bottom of previous layer
484        ktop=kbot
485        if     (ztop.ge.z_win(ktop+1)) then
486          ktop=ktop+1
487        endif
488       
489        zbot=z_wout(k+1)
490        zthk=zbot-ztop
491
492        if     (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then
493
494          kbot=ktop
495          do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
496            kbot=kbot+1
497          enddo
498          zbox=zbot
499          do k1= k+1,Nout
500            if     (z_wout(k1+1)-z_wout(k1).gt.dpthin) then
501              exit !thick layer
502            else
503              zbox=z_wout(k1+1)  !include thin adjacent layers
504              if     (zbox.eq.z_wout(Nout+1)) then
505                exit !at bottom
506              endif
507            endif
508          enddo
509          zthk=zbox-ztop
510
511          kbox=ktop
512          do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
513            kbox=kbox+1
514          enddo
515         
516          if     (ktop.eq.kbox) then
517
518
519            if     (z_wout(k)  .ne.z_win(kbox).or.z_wout(k+1).ne.z_win(kbox+1)     ) then
520
521              if     (hin(kbox).gt.dpthin) then
522                q001 = (zbox-z_win(kbox))/hin(kbox)
523                q002 = (ztop-z_win(kbox))/hin(kbox)
524                q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
525                q02=q01-1.+(q001+q002)
526                q0=1.-q01-q02
527              else
528                q0 = 1.0
529                q01 = 0.
530                q02 = 0.
531              endif
532          tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2) &
533       +q02*coeffremap(kbox,3)
534            else
535            tabout(k) = tabin(kbox)
536             
537            endif
538
539          else
540
541            if     (ktop.le.k .and. kbox.ge.k) then
542              ka = k
543            elseif (kbox-ktop.ge.3) then
544              ka = (kbox+ktop)/2
545            elseif (hin(ktop).ge.hin(kbox)) then
546              ka = ktop
547            else
548              ka = kbox
549            endif !choose ka
550
551            offset=coeffremap(ka,1)
552
553            qtop = z_win(ktop+1)-ztop !partial layer thickness
554            if     (hin(ktop).gt.dpthin) then
555              q=(ztop-z_win(ktop))/hin(ktop)
556              q01=q*(q-1.)
557              q02=q01+q
558              q0=1-q01-q02           
559            else
560              q0 = 1.
561              q01 = 0.
562              q02 = 0.
563            endif
564           
565            tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+ &
566       q02*coeffremap(ktop,3))-offset)*qtop
567
568            do k1= ktop+1,kbox-1
569              tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
570            enddo !k1
571
572           
573            qbot = zbox-z_win(kbox) !partial layer thickness
574            if     (hin(kbox).gt.dpthin) then
575              q=qbot/hin(kbox)
576              q01=(q-1.)**2
577              q02=q01-1.+q
578              q0=1-q01-q02                           
579            else
580              q0 = 1.0
581              q01 = 0.
582              q02 = 0.
583            endif
584           
585        tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+ &
586       q02*coeffremap(kbox,3))-offset)*qbot
587           
588            rpsum=1.0d0/zthk
589              tabout(k)=offset+tsum*rpsum
590             
591          endif !single or multiple layers
592        else
593        if (k==1) then
594        print *,'problem = ',zthk,z_wout(k+1),hout(1),hin(1),hin(2)
595        stop
596        endif
597         tabout(k) = tabout(k-1)
598         
599        endif !normal:thin layer
600      enddo !k
601           
602      return
603!---------------------------------------------------------------------------------------------------
604      end subroutine reconstructandremap_ps
605!===================================================================================================
606!
607
608end module module_interp
Note: See TracBrowser for help on using the repository browser.