New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
vremap.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/vremap.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
File size: 15.5 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, remap_linear
23
24   !! * Substitutions
25#  include "do_loop_substitute.h90"
26   !!----------------------------------------------------------------------
27   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
28   !! $Id: vremap 11573 2019-09-19 09:18:03Z jchanut $
29   !! Software governed by the CeCILL license (see ./LICENSE)
30   !!----------------------------------------------------------------------
31CONTAINS
32
33#if ! defined PPR_LIB
34   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var)     
35      !!----------------------------------------------------------------------
36      !!                    ***  ROUTINE  reconstructandremap ***
37      !!
38      !! ** Purpose :   Brief description of the routine
39      !!
40      !! ** Method  :   description of the methodoloy used to achieve the
41      !!                objectives of the routine. Be as clear as possible!
42      !!
43      !! ** Action  : - first action (share memory array/varible modified
44      !!                in this routine
45      !!              - second action .....
46      !!              - .....
47      !!
48      !! References :   Author et al., Short_name_review, Year
49      !!                Give references if exist otherwise suppress these lines
50      !!-----------------------------------------------------------------------
51      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
52      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
53      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
54      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   phin       ! Input thicknesses
55      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   phout      ! Output thicknesses
56      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
57      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Remapped data
58      !
59      INTEGER             :: jk, jn, k1, kbox, ktop, ka, kbot
60      REAL(wp), PARAMETER :: dpthin = 1.D-3, dsmll = 1.0D-8
61      REAL(wp)            :: q, q01, q02, q001, q002, q0
62      REAL(wp)            :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop
63      REAL(wp)            :: coeffremap(kjpk_in,3), zwork(kjpk_in,3), zwork2(kjpk_in+1,3)
64      REAL(wp)            :: z_win(1:kjpk_in+1), z_wout(1:kjpk_out+1)
65      !!-----------------------------------------------------------------------
66
67      z_win(1)=0._wp ; z_wout(1)= 0._wp
68      DO jk = 1, kjpk_in
69         z_win(jk+1)=z_win(jk)+phin(jk)
70      END DO
71     
72      DO jk = 1, kjpk_out
73         z_wout(jk+1)=z_wout(jk)+phout(jk)       
74      END DO       
75
76      DO jk = 2, kjpk_in
77         zwork(jk,1)=1._wp/(phin(jk-1)+phin(jk))
78      END DO
79       
80      DO jk = 2, kjpk_in-1
81         q0 = 1._wp / (phin(jk-1)+phin(jk)+phin(jk+1))
82         zwork(jk,2) = phin(jk-1) + 2._wp*phin(jk) + phin(jk+1)
83         zwork(jk,3) = q0
84      END DO       
85     
86      DO jn = 1, kn_var
87
88         DO jk = 2,kjpk_in
89            zwork2(jk,1) = zwork(jk,1)*(ptin(jk,jn)-ptin(jk-1,jn))
90         END DO
91       
92         coeffremap(:,1) = ptin(:,jn)
93 
94         DO jk = 2, kjpk_in-1
95            q001 = phin(jk)*zwork2(jk+1,1)
96            q002 = phin(jk)*zwork2(jk,1)       
97            IF (q001*q002 < 0._wp) then
98               q001 = 0._wp
99               q002 = 0._wp
100            ENDIF
101            q=zwork(jk,2)
102            q01=q*zwork2(jk+1,1)
103            q02=q*zwork2(jk,1)
104            IF (abs(q001) > abs(q02)) q001 = q02
105            IF (abs(q002) > abs(q01)) q002 = q01
106       
107            q=(q001-q002)*zwork(jk,3)
108            q001=q001-q*phin(jk+1)
109            q002=q002+q*phin(jk-1)
110       
111            coeffremap(jk,3)=coeffremap(jk,1)+q001
112            coeffremap(jk,2)=coeffremap(jk,1)-q002
113       
114            zwork2(jk,1)=(2._wp*q001-q002)**2
115            zwork2(jk,2)=(2._wp*q002-q001)**2
116         ENDDO
117       
118         DO jk = 1, kjpk_in
119            IF(jk.EQ.1 .OR. jk.EQ.kjpk_in .OR. phin(jk).LE.dpthin) THEN
120               coeffremap(jk,3) = coeffremap(jk,1)
121               coeffremap(jk,2) = coeffremap(jk,1)
122               zwork2(jk,1) = 0._wp
123               zwork2(jk,2) = 0._wp
124            ENDIF
125         END DO
126       
127         DO jk = 2, kjpk_in
128            q002 = max(zwork2(jk-1,2),dsmll)
129            q001 = max(zwork2(jk,1)  ,dsmll)
130            zwork2(jk,3) = (q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002)
131         END DO
132       
133         zwork2(1,3) = 2._wp*coeffremap(1,1)-zwork2(2,3)
134         zwork2(kjpk_in+1,3)=2._wp*coeffremap(kjpk_in,1)-zwork2(kjpk_in,3)
135 
136         DO jk = 1, kjpk_in
137            q01=zwork2(jk+1,3)-coeffremap(jk,1)
138            q02=coeffremap(jk,1)-zwork2(jk,3)
139            q001=2._wp*q01
140            q002=2._wp*q02
141            IF (q01*q02<0._wp) then
142               q01=0._wp
143               q02=0._wp
144            ELSEIF (abs(q01)>abs(q002)) then
145               q01=q002
146            ELSEIF (abs(q02)>abs(q001)) then
147               q02=q001
148            ENDIF
149            coeffremap(jk,2)=coeffremap(jk,1)-q02
150            coeffremap(jk,3)=coeffremap(jk,1)+q01
151         ENDDO
152
153         zbot=0._wp
154         kbot=1
155         DO jk=1,kjpk_out
156            ztop=zbot  !top is bottom of previous layer
157            ktop=kbot
158            IF     (ztop.GE.z_win(ktop+1)) then
159               ktop=ktop+1
160            ENDIF
161       
162            zbot=z_wout(jk+1)
163            zthk=zbot-ztop
164
165            IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(kjpk_out+1)) THEN
166 
167               kbot=ktop
168               DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.kjpk_in)
169                  kbot=kbot+1
170               ENDDO
171               zbox=zbot
172               DO k1= jk+1,kjpk_out
173                  IF     (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN
174                     exit !thick layer
175                  ELSE
176                     zbox=z_wout(k1+1)  !include thin adjacent layers
177                     IF(zbox.EQ.z_wout(kjpk_out+1)) THEN
178                        exit !at bottom
179                     ENDIF
180                  ENDIF
181               ENDDO
182               zthk=zbox-ztop
183
184               kbox=ktop
185               DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.kjpk_in)
186                  kbox=kbox+1
187               ENDDO
188         
189               IF(ktop.EQ.kbox) THEN
190                  IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN
191                     IF(phin(kbox).GT.dpthin) THEN
192                        q001 = (zbox-z_win(kbox))/phin(kbox)
193                        q002 = (ztop-z_win(kbox))/phin(kbox)
194                        q01=q001**2+q002**2+q001*q002+1._wp-2._wp*(q001+q002)
195                        q02=q01-1._wp+(q001+q002)
196                        q0=1._wp-q01-q02
197                     ELSE
198                        q0  = 1._wp
199                        q01 = 0._wp
200                        q02 = 0._wp
201                     ENDIF
202                     ptout(jk,jn)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3)
203                  ELSE
204                     ptout(jk,jn) = ptin(kbox,jn)
205                  ENDIF
206               ELSE
207                  IF(ktop.LE.jk .AND. kbox.GE.jk) THEN
208                     ka = jk
209                  ELSEIF (kbox-ktop.GE.3) THEN
210                     ka = (kbox+ktop)/2
211                  ELSEIF (phin(ktop).GE.phin(kbox)) THEN
212                     ka = ktop
213                  ELSE
214                     ka = kbox
215                  ENDIF !choose ka
216   
217                  offset=coeffremap(ka,1)
218   
219                  qtop = z_win(ktop+1)-ztop !partial layer thickness
220                  IF(phin(ktop).GT.dpthin) THEN
221                     q=(ztop-z_win(ktop))/phin(ktop)
222                     q01=q*(q-1._wp)
223                     q02=q01+q
224                     q0=1._wp-q01-q02           
225                  ELSE
226                     q0  = 1._wp
227                     q01 = 0._wp
228                     q02 = 0._wp
229                  ENDIF
230               
231                  tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop
232   
233                  DO k1= ktop+1,kbox-1
234                     tsum =tsum +(coeffremap(k1,1)-offset)*phin(k1)
235                  ENDDO !k1
236               
237                  qbot = zbox-z_win(kbox) !partial layer thickness
238                  IF(phin(kbox).GT.dpthin) THEN
239                     q=qbot/phin(kbox)
240                     q01=(q-1._wp)**2
241                     q02=q01-1._wp+q
242                     q0=1_wp-q01-q02                           
243                  ELSE
244                     q0  = 1._wp
245                     q01 = 0._wp
246                     q02 = 0._wp
247                  ENDIF
248             
249                  tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot
250               
251                  rpsum=1._wp / zthk
252                  ptout(jk,jn)=offset+tsum*rpsum
253                 
254               ENDIF !single or multiple layers
255            ELSE
256               IF (jk==1) THEN
257                  write(*,'(a7,i4,i4,3f12.5)')'problem = ',kjpk_in,kjpk_out,zthk,z_wout(jk+1),phout(1)
258               ENDIF
259               ptout(jk,jn) = ptout(jk-1,jn)
260             
261            ENDIF !normal:thin layer
262         ENDDO !jk
263
264      END DO ! loop over variables
265           
266   END SUBROUTINE reconstructandremap
267
268#else
269
270   SUBROUTINE reconstructandremap(ptin, phin, ptout, phout, kjpk_in, kjpk_out, kn_var)
271      !!----------------------------------------------------------------------
272      !!                    *** ROUTINE  reconstructandremap ***
273      !!
274      !! ** Purpose :   Conservative remapping of a vertical column
275      !!                from one set of layers to an other one.
276      !!
277      !! ** Method  :   Uses D. Engwirda Piecewise Polynomial Reconstruction library.
278      !!                https://github.com/dengwirda/PPR
279      !!               
280      !!
281      !! References :   Engwirda, Darren & Kelley, Maxwell. (2015). A WENO-type
282      !!                slope-limiter for a family of piecewise polynomial methods.
283      !!                https://arxiv.org/abs/1606.08188
284      !!-----------------------------------------------------------------------
285      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
286      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
287      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
288      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   phin       ! Input thicknesses
289      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   phout      ! Output thicknesses
290      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
291      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Remapped data
292      !
293      INTEGER, PARAMETER :: ndof = 1
294      INTEGER  :: jk, jn
295      REAL(wp) ::  zwin(kjpk_in+1) ,  ztin(ndof, kn_var, kjpk_in)
296      REAL(wp) :: zwout(kjpk_out+1), ztout(ndof, kn_var, kjpk_out)
297      TYPE(rmap_work) :: work
298      TYPE(rmap_opts) :: opts
299      TYPE(rcon_ends) :: bc_l(kn_var)
300      TYPE(rcon_ends) :: bc_r(kn_var)
301      !!--------------------------------------------------------------------
302     
303      ! Set interfaces and input data:
304      zwin(1) = 0._wp
305      DO jk = 2, kjpk_in + 1
306         zwin(jk) = zwin(jk-1) + phin(jk-1) 
307      END DO
308     
309      DO jn = 1, kn_var 
310         DO jk = 1, kjpk_in
311            ztin(ndof, jn, jk) =  ptin(jk, jn)
312         END DO
313      END DO
314
315      zwout(1) = 0._wp
316      DO jk = 2, kjpk_out + 1
317         zwout(jk) = zwout(jk-1) + phout(jk-1) 
318      END DO
319
320      ! specify methods
321!      opts%edge_meth = p1e_method     ! 1st-order edge interp.
322!      opts%cell_meth = plm_method     ! PLM method in cells
323      opts%edge_meth = p3e_method     ! 3rd-order edge interp.
324      opts%cell_meth = ppm_method     ! PPM method in cells   
325!      opts%edge_meth = p5e_method     ! 5th-order edge interp.
326!      opts%cell_meth = pqm_method     ! PQM method in cells
327
328      ! limiter
329!      opts%cell_lims = null_limit     ! no lim.
330      opts%cell_lims = mono_limit     ! monotone limiter   
331 
332      ! set boundary conditions
333      bc_l%bcopt = bcon_loose         ! "loose" = extrapolate
334      bc_r%bcopt = bcon_loose
335!      bc_l%bcopt = bcon_slope       
336!      bc_r%bcopt = bcon_slope
337
338      ! init. method workspace
339      CALL work%init(kjpk_in+1, kn_var, opts)
340
341      ! remap
342      CALL rmap1d(kjpk_in+1, kjpk_out+1, kn_var, ndof, &
343      &           zwin, zwout, ztin, ztout,            &
344      &           bc_l, bc_r, work, opts)
345
346      ! clear method workspace
347      CALL work%free()
348
349      DO jn = 1, kn_var 
350         DO jk = 1, kjpk_out
351            ptout(jk, jn) = ztout(1, jn, jk)
352         END DO
353      END DO
354           
355   END SUBROUTINE reconstructandremap
356#endif
357
358   SUBROUTINE remap_linear(ptin, pzin, ptout, pzout, kjpk_in, kjpk_out, kn_var)
359      !!----------------------------------------------------------------------
360      !!                    *** ROUTINE  remap_linear ***
361      !!
362      !! ** Purpose :   Linear interpolation based on input/ouputs depths
363      !!
364      !!-----------------------------------------------------------------------
365      INTEGER , INTENT(in   )                      ::   kjpk_in    ! Number of input levels
366      INTEGER , INTENT(in   )                      ::   kjpk_out   ! Number of output levels
367      INTEGER , INTENT(in   )                      ::   kn_var     ! Number of variables
368      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in)  ::   pzin       ! Input depths
369      REAL(wp), INTENT(in   ), DIMENSION(kjpk_out) ::   pzout      ! Output depths
370      REAL(wp), INTENT(in   ), DIMENSION(kjpk_in , kn_var) ::   ptin       ! Input data
371      REAL(wp), INTENT(inout), DIMENSION(kjpk_out, kn_var) ::   ptout      ! Interpolated data
372      !
373      INTEGER  :: jkin, jkout, jn
374      !!--------------------------------------------------------------------
375      !     
376      DO jkout = 1, kjpk_out !  Loop over destination grid
377         !
378         IF     ( pzout(jkout) <=  pzin(  1    ) ) THEN  ; ptout(jkout,1:kn_var) = ptin(    1    ,1:kn_var)         
379         ELSEIF ( pzout(jkout) >= pzin(kjpk_in) ) THEN   ; ptout(jkout,1:kn_var) = ptin( kjpk_in ,1:kn_var)
380         ELSEIF ( ( pzout(jkout) > pzin(1) ).AND.( pzout(jkout) < pzin(kjpk_in) )) THEN
381            DO jkin = 1, kjpk_in - 1 !  Loop over source grid
382               IF ( pzout(jkout) < pzin(jkin+1) ) THEN
383                  DO jn = 1, kn_var
384                     ptout(jkout,jn) =  ptin(jkin,jn) + &
385                                     & (pzout(jkout) - pzin(jkin)) / (pzin(jkin+1)    - pzin(jkin)) &
386                                     &                             * (ptin(jkin+1,jn) - ptin(jkin,jn))
387                  END DO 
388                  EXIT
389               ENDIF 
390            END DO
391         ENDIF
392         !
393      END DO
394
395   END SUBROUTINE remap_linear
396
397   !!======================================================================
398!$AGRIF_END_DO_NOT_TREAT
399END MODULE vremap
Note: See TracBrowser for help on using the repository browser.