source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/vremap.F90 @ 11590

Last change on this file since 11590 was 11590, checked in by jchanut, 16 months ago

#2222: 1) create remapping module (vremap) and integration of D. Engwirda piecewise polynomial recontruction package (PPR_LIB cpp key). 2) Various bug corrections with key_vertical activated.

File size: 11.4 KB
Line 
1#undef PPR_LIB      /* USE PPR library */
2MODULE vremap
3!$AGRIF_DO_NOT_TREAT
4   !!======================================================================
5   !!                       ***  MODULE  vremap  ***
6   !! Ocean physics:  Vertical remapping routines
7   !!
8   !!======================================================================
9   !! History : 4.0  !  2019-09  (Jérôme Chanut)  Original code
10   !!----------------------------------------------------------------------
11   !!----------------------------------------------------------------------
12   !!
13   !!----------------------------------------------------------------------
14   USE par_oce
15#if defined PPR_LIB
16   USE ppr_1d   ! D. Engwirda piecewise polynomial reconstruction library
17#endif
18
19   IMPLICIT NONE
20   PRIVATE
21
22   PUBLIC   reconstructandremap
23
24   !!----------------------------------------------------------------------
25   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
26   !! $Id: vremap 11573 2019-09-19 09:18:03Z jchanut $
27   !! Software governed by the CeCILL license (see ./LICENSE)
28   !!----------------------------------------------------------------------
29CONTAINS
30
31#if ! defined PPR_LIB
32   SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout)     
33      !!----------------------------------------------------------------------
34      !!                    ***  ROUTINE  reconstructandremap ***
35      !!
36      !! ** Purpose :   Brief description of the routine
37      !!
38      !! ** Method  :   description of the methodoloy used to achieve the
39      !!                objectives of the routine. Be as clear as possible!
40      !!
41      !! ** Action  : - first action (share memory array/varible modified
42      !!                in this routine
43      !!              - second action .....
44      !!              - .....
45      !!
46      !! References :   Author et al., Short_name_review, Year
47      !!                Give references if exist otherwise suppress these lines
48      !!-----------------------------------------------------------------------
49      INTEGER N, Nout
50      REAL(wp) tabin(N), tabout(Nout)
51      REAL(wp) hin(N), hout(Nout)
52      REAL(wp) coeffremap(N,3),zwork(N,3)
53      REAL(wp) zwork2(N+1,3)
54      INTEGER jk
55      DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8 
56      REAL(wp) q,q01,q02,q001,q002,q0
57      REAL(wp) z_win(1:N+1), z_wout(1:Nout+1)
58      REAL(wp),PARAMETER :: dpthin = 1.D-3
59      INTEGER :: k1, kbox, ktop, ka, kbot
60      REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
61      !!-----------------------------------------------------------------------
62
63      z_win(1)=0.; z_wout(1)= 0.
64      DO jk=1,N
65         z_win(jk+1)=z_win(jk)+hin(jk)
66      ENDDO 
67     
68      DO jk=1,Nout
69         z_wout(jk+1)=z_wout(jk)+hout(jk)       
70      ENDDO       
71
72      DO jk=2,N
73         zwork(jk,1)=1./(hin(jk-1)+hin(jk))
74      ENDDO
75       
76      DO jk=2,N-1
77         q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1))
78         zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1)
79         zwork(jk,3)=q0
80      ENDDO       
81     
82      DO jk= 2,N
83         zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1))
84      ENDDO
85       
86      coeffremap(:,1) = tabin(:)
87 
88      DO jk=2,N-1
89         q001 = hin(jk)*zwork2(jk+1,1)
90         q002 = hin(jk)*zwork2(jk,1)       
91         IF (q001*q002 < 0) then
92            q001 = 0.
93            q002 = 0.
94         ENDIF
95         q=zwork(jk,2)
96         q01=q*zwork2(jk+1,1)
97         q02=q*zwork2(jk,1)
98         IF (abs(q001) > abs(q02)) q001 = q02
99         IF (abs(q002) > abs(q01)) q002 = q01
100       
101         q=(q001-q002)*zwork(jk,3)
102         q001=q001-q*hin(jk+1)
103         q002=q002+q*hin(jk-1)
104       
105         coeffremap(jk,3)=coeffremap(jk,1)+q001
106         coeffremap(jk,2)=coeffremap(jk,1)-q002
107       
108         zwork2(jk,1)=(2.*q001-q002)**2
109         zwork2(jk,2)=(2.*q002-q001)**2
110      ENDDO
111       
112      DO jk=1,N
113         IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN
114            coeffremap(jk,3) = coeffremap(jk,1)
115            coeffremap(jk,2) = coeffremap(jk,1)
116            zwork2(jk,1) = 0.
117            zwork2(jk,2) = 0.
118         ENDIF
119      ENDDO
120       
121      DO jk=2,N
122         q002=max(zwork2(jk-1,2),dsmll)
123         q001=max(zwork2(jk,1),dsmll)
124         zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)
125      ENDDO
126       
127      zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3)
128      zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3)
129 
130      DO jk=1,N
131         q01=zwork2(jk+1,3)-coeffremap(jk,1)
132         q02=coeffremap(jk,1)-zwork2(jk,3)
133         q001=2.*q01
134         q002=2.*q02
135         IF (q01*q02<0) then
136            q01=0.
137            q02=0.
138         ELSEIF (abs(q01)>abs(q002)) then
139            q01=q002
140         ELSEIF (abs(q02)>abs(q001)) then
141            q02=q001
142         ENDIF
143         coeffremap(jk,2)=coeffremap(jk,1)-q02
144         coeffremap(jk,3)=coeffremap(jk,1)+q01
145      ENDDO
146
147      zbot=0.0
148      kbot=1
149      DO jk=1,Nout
150         ztop=zbot  !top is bottom of previous layer
151         ktop=kbot
152         IF     (ztop.GE.z_win(ktop+1)) then
153            ktop=ktop+1
154         ENDIF
155       
156         zbot=z_wout(jk+1)
157         zthk=zbot-ztop
158
159         IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN
160
161            kbot=ktop
162            DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N)
163               kbot=kbot+1
164            ENDDO
165            zbox=zbot
166            DO k1= jk+1,Nout
167               IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
168                  exit !thick layer
169               ELSE
170                  zbox=z_wout(k1+1)  !include thin adjacent layers
171                  IF(zbox.EQ.z_wout(Nout+1)) THEN
172                     exit !at bottom
173                  ENDIF
174               ENDIF
175            ENDDO
176            zthk=zbox-ztop
177
178            kbox=ktop
179            DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N)
180               kbox=kbox+1
181            ENDDO
182         
183            IF(ktop.EQ.kbox) THEN
184               IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
185                  IF(hin(kbox).GT.dpthin) THEN
186                     q001 = (zbox-z_win(kbox))/hin(kbox)
187                     q002 = (ztop-z_win(kbox))/hin(kbox)
188                     q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002)
189                     q02=q01-1.+(q001+q002)
190                     q0=1.-q01-q02
191                  ELSE
192                     q0 = 1.0
193                     q01 = 0.
194                     q02 = 0.
195                  ENDIF
196                  tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
197               ELSE
198                  tabout(jk) = tabin(kbox)
199               ENDIF
200            ELSE
201               IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
202                  ka = jk
203               ELSEIF (kbox-ktop.GE.3) THEN
204                  ka = (kbox+ktop)/2
205               ELSEIF (hin(ktop).GE.hin(kbox)) THEN
206                  ka = ktop
207               ELSE
208                  ka = kbox
209               ENDIF !choose ka
210   
211               offset=coeffremap(ka,1)
212   
213               qtop = z_win(ktop+1)-ztop !partial layer thickness
214               IF(hin(ktop).GT.dpthin) THEN
215                  q=(ztop-z_win(ktop))/hin(ktop)
216                  q01=q*(q-1.)
217                  q02=q01+q
218                  q0=1-q01-q02           
219               ELSE
220                  q0 = 1.
221                  q01 = 0.
222                  q02 = 0.
223               ENDIF
224               
225               tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
226   
227               DO k1= ktop+1,kbox-1
228                  tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1)
229               ENDDO !k1
230               
231               qbot = zbox-z_win(kbox) !partial layer thickness
232               IF(hin(kbox).GT.dpthin) THEN
233                  q=qbot/hin(kbox)
234                  q01=(q-1.)**2
235                  q02=q01-1.+q
236                  q0=1-q01-q02                           
237               ELSE
238                  q0 = 1.0
239                  q01 = 0.
240                  q02 = 0.
241               ENDIF
242             
243               tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
244               
245               rpsum=1.0d0/zthk
246               tabout(jk)=offset+tsum*rpsum
247                 
248            ENDIF !single or multiple layers
249         ELSE
250            IF (jk==1) THEN
251               write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1)
252            ENDIF
253            tabout(jk) = tabout(jk-1)
254             
255         ENDIF !normal:thin layer
256      ENDDO !jk
257           
258   END SUBROUTINE reconstructandremap
259
260#else
261
262   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, jkin, jkout)
263      !!----------------------------------------------------------------------
264      !!                    ***  ROUTINE  reconstructandremap ***
265      !!
266      !! ** Purpose :   Brief description of the routine
267      !!
268      !! ** Method  :   description of the methodoloy used to achieve the
269      !!                objectives of the routine. Be as clear as possible!
270      !!
271      !! ** Action  : - first action (share memory array/varible modified
272      !!                in this routine
273      !!              - second action .....
274      !!              - .....
275      !!
276      !! References :   Author et al., Short_name_review, Year
277      !!                Give references if exist otherwise suppress these lines
278      !!-----------------------------------------------------------------------
279      INTEGER , INTENT(in   )                     ::   jkin, jkout       !
280      REAL(wp), INTENT(in   ), DIMENSION(jkin)    ::   phin, ptin        !
281      REAL(wp), INTENT(in   ), DIMENSION(jkout)   ::   phout             !
282      REAL(wp), INTENT(inout), DIMENSION(jkout)   ::   ptout             !
283      !
284      INTEGER, PARAMETER :: ndof = 1, nvar = 1
285      INTEGER  :: jk
286      REAL(wp) ::  zwin(jkin+1),   ztin(ndof, nvar, jkin)
287      REAL(wp) :: zwout(jkout+1), ztout(ndof, nvar, jkout)
288      TYPE(rmap_work) :: work
289      TYPE(rmap_opts) :: opts
290      TYPE(rcon_ends) :: bc_l(nvar)
291      TYPE(rcon_ends) :: bc_r(nvar)
292      !!--------------------------------------------------------------------
293     
294      ! Set interfaces and input data:
295      zwin(1) = 0._wp
296      DO jk = 2, jkin + 1
297         zwin(jk) = zwin(jk-1) + phin(jk-1) 
298         ztin(ndof, nvar,jk-1) =  ptin(jk-1)
299      END DO
300
301      zwout(1) = 0._wp
302      DO jk = 2, jkout + 1
303         zwout(jk) = zwout(jk-1) + phout(jk-1) 
304      END DO
305
306      ! specify methods
307!      opts%edge_meth = p3e_method     ! 3rd-order edge interp.
308!      opts%cell_meth = ppm_method     ! PPM method in cells
309!      opts%cell_lims = null_limit     ! no lim.     
310      opts%edge_meth = p5e_method     ! 5th-order edge interp.
311      opts%cell_meth = pqm_method     ! PPM method in cells
312      opts%cell_lims = mono_limit     ! monotone limiter   
313 
314      ! set boundary conditions
315      bc_l%bcopt = bcon_loose         ! "loose" = extrapolate
316      bc_r%bcopt = bcon_loose
317!      bc_l%bcopt = bcon_slope       
318!      bc_r%bcopt = bcon_slope
319
320      ! init. method workspace
321      CALL work%init(jkin+1, nvar, opts)
322
323      ! remap
324      CALL rmap1d(jkin+1, jkout+1, nvar, ndof, &
325      &           zwin, zwout, ztin, ztout,    &
326      &           bc_l, bc_r, work, opts)
327
328      ! clear method workspace
329      CALL work%free()
330
331      DO jk =1, jkout
332         ptout(jk) = ztout(1,1,jk)
333      END DO
334           
335   END SUBROUTINE reconstructandremap
336#endif
337
338   !!======================================================================
339!$AGRIF_END_DO_NOT_TREAT
340END MODULE vremap
Note: See TracBrowser for help on using the repository browser.