source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF/isftbl.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 8 months 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: 12.3 KB
Line 
1MODULE isftbl
2   !!======================================================================
3   !!                       ***  MODULE  isftbl  ***
4   !! isftbl module :  compute properties of top boundary layer
5   !!======================================================================
6   !! History :  4.1  !  2019-09  (P. Mathiot) original code
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   isftbl       : routine to compute :
11   !!                  - geometry of the ice shelf tbl (isf_tbl_lvl, isftbl_ktop, isftbl_kbot)
12   !!                    (top and bottom level, thickness and fraction of deepest level affected)
13   !!                  - tbl averaged properties (isf_tbl, isf_tbl_avg)
14   !!----------------------------------------------------------------------
15
16   USE isf_oce ! ice shelf variables
17
18   USE dom_oce ! vertical scale factor and depth
19
20   IMPLICIT NONE
21
22   PRIVATE
23
24   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot
25   !! * Substitutions
26#  include "do_loop_substitute.h90"
27
28CONTAINS
29
30   SUBROUTINE isf_tbl( Kmm, pvarin, pvarout, cd_ptin, ktop, phtbl, kbot, pfrac )
31      !!--------------------------------------------------------------------
32      !!                  ***  SUBROUTINE isf_tbl  ***
33      !!
34      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
35      !!
36      !! ** Method : Average properties over a specific thickness
37      !!
38      !! ** Reference : inspired from : Losch, Modeling ice shelf cavities in a z coordinate ocean general circulation model
39      !!                https://doi.org/10.1029/2007JC004368 , 2008
40      !!
41      !!-------------------------- OUT -------------------------------------
42      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) :: pvarout ! 2d average of pvarin
43      !!-------------------------- IN  -------------------------------------
44      INTEGER                               , INTENT(in   ) :: Kmm           ! ocean time level index
45      CHARACTER(len=1)                      , INTENT(in   ) :: cd_ptin       ! point of variable in/out
46      REAL(wp), DIMENSION(jpi,jpj,jpk)      , INTENT(in   ) :: pvarin        ! 3d variable to average over the tbl
47      INTEGER,  DIMENSION(jpi,jpj)          , INTENT(in   ) :: ktop          ! top level
48      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in   ) :: phtbl         ! tbl thickness
49      !!-------------------------- IN OPTIONAL -----------------------------
50      INTEGER,  DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: kbot          ! bottom level
51      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(in   ) :: pfrac         ! fraction of bottom cell affected by tbl
52      !!--------------------------------------------------------------------
53      INTEGER ::   ji, jj                     ! loop index
54      INTEGER , DIMENSION(jpi,jpj) :: ikbot   ! bottom level of the tbl
55      REAL(wp), DIMENSION(jpi,jpj) :: zvarout ! 2d average of pvarin
56      REAL(wp), DIMENSION(jpi,jpj) :: zhtbl   ! thickness of the tbl
57      REAL(wp), DIMENSION(jpi,jpj) :: zfrac   ! thickness of the tbl
58      !!--------------------------------------------------------------------
59      !
60      SELECT CASE ( cd_ptin )
61      CASE ( 'U' )
62         !
63         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
64         zhtbl = phtbl
65         !
66         ! compute tbl lvl and thickness
67         CALL isf_tbl_lvl( hu(:,:,Kmm), e3u(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
68         !
69         ! compute tbl property at U point
70         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u(:,:,:,Kmm), pvarin, zvarout )
71         !
72         ! compute tbl property at T point
73         pvarout(1,:) = 0._wp
74         DO_2D_11_01
75            pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
76         END_2D
77         ! lbclnk not needed as a final communication is done after the computation of fwf
78         !
79      CASE ( 'V' )
80         !
81         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
82         zhtbl = phtbl
83         !
84         ! compute tbl lvl and thickness
85         CALL isf_tbl_lvl( hv(:,:,Kmm), e3v(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
86         !
87         ! compute tbl property at V point
88         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v(:,:,:,Kmm), pvarin, zvarout )
89         !
90         ! pvarout is an averaging of wet point
91         pvarout(:,1) = 0._wp
92         DO_2D_01_11
93            pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
94         END_2D
95         ! lbclnk not needed as a final communication is done after the computation of fwf
96         !
97      CASE ( 'T' )
98         !
99         ! compute tbl property at T point
100         CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, e3t(:,:,:,Kmm), pvarin, pvarout )
101         !
102      END SELECT
103      !
104   END SUBROUTINE isf_tbl
105
106   SUBROUTINE isf_tbl_avg( ktop, kbot, phtbl, pfrac, pe3, pvarin, pvarout )
107      !!--------------------------------------------------------------------
108      !!                  ***  ROUTINE isf_tbl_avg  ***
109      !!
110      !! ** Purpose : compute mean property in the boundary layer
111      !!
112      !! ** Method  : Depth average is made between the top level ktop and the bottom level kbot
113      !!              over a thickness phtbl. The bottom level is partially counted (pfrac).
114      !!
115      !!-------------------------- OUT -------------------------------------
116      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pvarout      ! tbl property averaged over phtbl between level ktop and kbot
117      !!-------------------------- IN  -------------------------------------
118      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop, kbot   ! top and bottom level of the top boundary layer
119      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl, pfrac ! fraction of bottom level to be affected by the tbl
120      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3          ! vertical scale factor
121      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pvarin       ! tbl property to average between ktop, kbot over phtbl
122      !!--------------------------------------------------------------------
123      INTEGER  :: ji,jj,jk                    ! loop indices
124      INTEGER  :: ikt, ikb                    ! top and bottom levels
125      !!--------------------------------------------------------------------
126      !
127      ! compute tbl top.bottom level and thickness
128      DO_2D_11_11
129         !
130         ! tbl top/bottom indices initialisation
131         ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
132         !
133         ! level fully include in the ice shelf boundary layer
134         pvarout(ji,jj) = SUM( pvarin(ji,jj,ikt:ikb-1) * pe3(ji,jj,ikt:ikb-1) ) / phtbl(ji,jj)
135         !
136         ! level partially include in ice shelf boundary layer
137         pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * pe3(ji,jj,ikb) / phtbl(ji,jj) * pfrac(ji,jj)
138         !
139      END_2D
140
141   END SUBROUTINE isf_tbl_avg
142
143   SUBROUTINE isf_tbl_lvl( phw, pe3, ktop, kbot, phtbl, pfrac )
144      !!--------------------------------------------------------------------
145      !!                  ***  ROUTINE isf_tbl_lvl  ***
146      !!
147      !! ** Purpose : - compute bottom level off the top boundary layer
148      !!              - thickness of the top boundary layer
149      !!              - fraction of the bottom level affected by the tbl
150      !!
151      !!-------------------------- OUT --------------------------------------
152      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
153      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pfrac  ! fraction of bottom level in the tbl
154      !!-------------------------- IN  --------------------------------------
155      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
156      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phw    ! water column thickness
157      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
158      !!-------------------------- INOUT ------------------------------------
159      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(inout) :: phtbl  ! top boundary layer thickness
160      !!---------------------------------------------------------------------
161      INTEGER :: ji,jj,jk
162      INTEGER :: ikt, ikb
163      !!---------------------------------------------------------------------
164      !
165      ! get htbl
166      DO_2D_11_11
167         !
168         ! tbl top/bottom indices initialisation
169         ikt = ktop(ji,jj)
170         !
171         ! limit the tbl to water thickness.
172         phtbl(ji,jj) = MIN( phtbl(ji,jj), phw(ji,jj) )
173         !
174         ! thickness of boundary layer must be at least the top level thickness
175         phtbl(ji,jj) = MAX( phtbl(ji,jj), pe3(ji,jj,ikt) )
176         !
177      END_2D
178      !
179      ! get ktbl
180      CALL isf_tbl_kbot(ktop, phtbl, pe3, kbot)
181      !
182      ! get pfrac
183      DO_2D_11_11
184         !
185         ! tbl top/bottom indices initialisation
186         ikt = ktop(ji,jj) ; ikb = kbot(ji,jj)
187         !
188         ! proportion of the bottom cell included in ice shelf boundary layer
189         pfrac(ji,jj) = ( phtbl(ji,jj) - SUM( pe3(ji,jj,ikt:ikb-1) ) ) / pe3(ji,jj,ikb)
190         !
191      END_2D
192      !
193   END SUBROUTINE isf_tbl_lvl
194   !
195   SUBROUTINE isf_tbl_kbot(ktop, phtbl, pe3, kbot)
196      !!--------------------------------------------------------------------
197      !!                  ***  ROUTINE isf_tbl_bot  ***
198      !!
199      !! ** Purpose : compute bottom level of the isf top boundary layer
200      !!
201      !!-------------------------- OUT -------------------------------------
202      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
203      !!-------------------------- IN  -------------------------------------
204      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) :: phtbl  ! top boundary layer thickness
205      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(in   ) :: ktop   ! top level of the top boundary layer
206      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) :: pe3    ! vertical scale factor
207      !!--------------------------------------------------------------------
208      INTEGER :: ji, jj
209      INTEGER :: ikt, ikb
210      !!--------------------------------------------------------------------
211      !
212      ! phtbl need to be bounded by water column thickness before
213      ! test: if htbl = water column thickness, should return mbathy
214      ! test: if htbl = 0 should return ktop (phtbl cap to e3t(ji,jj,1))
215      !
216      ! get ktbl
217      DO_2D_11_11
218         !
219         ! determine the deepest level influenced by the boundary layer
220         ikt = ktop(ji,jj)
221         ikb = ikt
222         DO WHILE ( SUM(pe3(ji,jj,ikt:ikb-1)) < phtbl(ji,jj ) ) ;  ikb = ikb + 1 ;  END DO
223         kbot(ji,jj) = ikb - 1
224         !
225      END_2D
226      !
227   END SUBROUTINE isf_tbl_kbot
228      !
229   SUBROUTINE isf_tbl_ktop(pdep, ktop)
230      !!--------------------------------------------------------------------
231      !!                  ***  ROUTINE isf_tbl_top  ***
232      !!
233      !! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
234      !!
235      !!-------------------------- OUT -------------------------------------
236      INTEGER,  DIMENSION(jpi,jpj), INTENT(  out) :: ktop        ! top level affected by the ice shelf parametrisation
237      !!-------------------------- IN  -------------------------------------
238      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdep        ! top depth of the parametrisation influence
239      !!--------------------------------------------------------------------
240      INTEGER :: ji,jj
241      INTEGER :: ikt
242      !!--------------------------------------------------------------------
243      !
244      ! if we need to recompute the top level at every time stepcompute top level (z*, z~)
245      ! in case of weak ht variation we can assume the top level of htbl to be constant
246      ! => only done using gdepw_0
247      ! be sure pdep is already correctly bounded
248      ! test: this routine run on isfdraft should return mikt
249      ! test: this routine run with pdep = 0 should return 1
250      !
251      DO_2D_11_11
252         ! comput ktop
253         ikt = 2
254         DO WHILE ( gdepw_0(ji,jj,ikt) <= pdep(ji,jj ) ) ;  ikt = ikt + 1 ;  END DO
255         ktop(ji,jj) = ikt - 1
256         !
257         ! update pdep
258         pdep(ji,jj) = gdepw_0(ji,jj,ktop(ji,jj))
259      END_2D
260      !
261   END SUBROUTINE isf_tbl_ktop
262
263END MODULE isftbl
Note: See TracBrowser for help on using the repository browser.