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.
isftbl.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/ISF – NEMO

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

Last change on this file since 12808 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: 12.3 KB
RevLine 
[11395]1MODULE isftbl
[11403]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   !!----------------------------------------------------------------------
[11395]8
[11403]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   !!----------------------------------------------------------------------
[11395]15
[12077]16   USE isf_oce ! ice shelf variables
[11852]17
[11876]18   USE dom_oce ! vertical scale factor and depth
[11395]19
[11403]20   IMPLICIT NONE
[11395]21
[11403]22   PRIVATE
[11395]23
[11541]24   PUBLIC isf_tbl, isf_tbl_avg, isf_tbl_lvl, isf_tbl_ktop, isf_tbl_kbot
[12340]25   !! * Substitutions
26#  include "do_loop_substitute.h90"
[11403]27
[11395]28CONTAINS
29
[12068]30   SUBROUTINE isf_tbl( Kmm, pvarin, pvarout, cd_ptin, ktop, phtbl, kbot, pfrac )
[11403]31      !!--------------------------------------------------------------------
32      !!                  ***  SUBROUTINE isf_tbl  ***
[11395]33      !!
34      !! ** Purpose : compute mean T/S/U/V in the boundary layer at T- point
35      !!
[11541]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      !!
[11395]41      !!-------------------------- OUT -------------------------------------
[11403]42      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(  out) :: pvarout ! 2d average of pvarin
[11395]43      !!-------------------------- IN  -------------------------------------
[12068]44      INTEGER                               , INTENT(in   ) :: Kmm           ! ocean time level index
[11403]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
[11495]47      INTEGER,  DIMENSION(jpi,jpj)          , INTENT(in   ) :: ktop          ! top level
48      REAL(wp), DIMENSION(jpi,jpj)          , INTENT(in   ) :: phtbl         ! tbl thickness
[11403]49      !!-------------------------- IN OPTIONAL -----------------------------
[11495]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
[11395]52      !!--------------------------------------------------------------------
[11844]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
[11403]58      !!--------------------------------------------------------------------
[11395]59      !
60      SELECT CASE ( cd_ptin )
61      CASE ( 'U' )
62         !
[11495]63         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
64         zhtbl = phtbl
[11395]65         !
[11495]66         ! compute tbl lvl and thickness
[12068]67         CALL isf_tbl_lvl( hu(:,:,Kmm), e3u(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
[11495]68         !
[11395]69         ! compute tbl property at U point
[12068]70         CALL isf_tbl_avg( miku, ikbot, zhtbl, zfrac, e3u(:,:,:,Kmm), pvarin, zvarout )
[11395]71         !
72         ! compute tbl property at T point
[11844]73         pvarout(1,:) = 0._wp
[12340]74         DO_2D_11_01
75            pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji-1,jj))
76         END_2D
[11844]77         ! lbclnk not needed as a final communication is done after the computation of fwf
[11395]78         !
79      CASE ( 'V' )
80         !
[11495]81         ! copy phtbl (phtbl is INTENT in as we don't want to change it)
82         zhtbl = phtbl
[11395]83         !
[11495]84         ! compute tbl lvl and thickness
[12068]85         CALL isf_tbl_lvl( hv(:,:,Kmm), e3v(:,:,:,Kmm), ktop, ikbot, zhtbl, zfrac )
[11495]86         !
[11395]87         ! compute tbl property at V point
[12068]88         CALL isf_tbl_avg( mikv, ikbot, zhtbl, zfrac, e3v(:,:,:,Kmm), pvarin, zvarout )
[11395]89         !
90         ! pvarout is an averaging of wet point
[11844]91         pvarout(:,1) = 0._wp
[12340]92         DO_2D_01_11
93            pvarout(ji,jj) = 0.5_wp * (zvarout(ji,jj) + zvarout(ji,jj-1))
94         END_2D
[11844]95         ! lbclnk not needed as a final communication is done after the computation of fwf
[11395]96         !
97      CASE ( 'T' )
98         !
99         ! compute tbl property at T point
[12068]100         CALL isf_tbl_avg( ktop, kbot, phtbl, pfrac, e3t(:,:,:,Kmm), pvarin, pvarout )
[11395]101         !
102      END SELECT
103      !
104   END SUBROUTINE isf_tbl
105
106   SUBROUTINE isf_tbl_avg( ktop, kbot, phtbl, pfrac, pe3, pvarin, pvarout )
[11403]107      !!--------------------------------------------------------------------
[11541]108      !!                  ***  ROUTINE isf_tbl_avg  ***
[11395]109      !!
110      !! ** Purpose : compute mean property in the boundary layer
111      !!
[11403]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      !!
[11395]115      !!-------------------------- OUT -------------------------------------
[11403]116      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pvarout      ! tbl property averaged over phtbl between level ktop and kbot
[11395]117      !!-------------------------- IN  -------------------------------------
[11403]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      !!--------------------------------------------------------------------
[11395]123      INTEGER  :: ji,jj,jk                    ! loop indices
124      INTEGER  :: ikt, ikb                    ! top and bottom levels
[11403]125      !!--------------------------------------------------------------------
[11395]126      !
127      ! compute tbl top.bottom level and thickness
[12340]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
[11395]140
141   END SUBROUTINE isf_tbl_avg
142
143   SUBROUTINE isf_tbl_lvl( phw, pe3, ktop, kbot, phtbl, pfrac )
[11403]144      !!--------------------------------------------------------------------
[11395]145      !!                  ***  ROUTINE isf_tbl_lvl  ***
146      !!
[11541]147      !! ** Purpose : - compute bottom level off the top boundary layer
[11395]148      !!              - thickness of the top boundary layer
[11541]149      !!              - fraction of the bottom level affected by the tbl
[11395]150      !!
[11495]151      !!-------------------------- OUT --------------------------------------
[11403]152      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
[11495]153      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(  out) :: pfrac  ! fraction of bottom level in the tbl
154      !!-------------------------- IN  --------------------------------------
[11403]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
[11495]158      !!-------------------------- INOUT ------------------------------------
159      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(inout) :: phtbl  ! top boundary layer thickness
[11395]160      !!---------------------------------------------------------------------
161      INTEGER :: ji,jj,jk
162      INTEGER :: ikt, ikb
163      !!---------------------------------------------------------------------
164      !
165      ! get htbl
[12340]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
[11395]178      !
179      ! get ktbl
[11541]180      CALL isf_tbl_kbot(ktop, phtbl, pe3, kbot)
[11395]181      !
182      ! get pfrac
[12340]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
[11395]192      !
193   END SUBROUTINE isf_tbl_lvl
194   !
[11541]195   SUBROUTINE isf_tbl_kbot(ktop, phtbl, pe3, kbot)
[11403]196      !!--------------------------------------------------------------------
[11541]197      !!                  ***  ROUTINE isf_tbl_bot  ***
[11403]198      !!
199      !! ** Purpose : compute bottom level of the isf top boundary layer
200      !!
[11395]201      !!-------------------------- OUT -------------------------------------
[11403]202      INTEGER,  DIMENSION(jpi,jpj)    , INTENT(  out) :: kbot   ! bottom level of the top boundary layer
[11395]203      !!-------------------------- IN  -------------------------------------
[11403]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      !!--------------------------------------------------------------------
[11395]208      INTEGER :: ji, jj
209      INTEGER :: ikt, ikb
[11403]210      !!--------------------------------------------------------------------
[11395]211      !
212      ! phtbl need to be bounded by water column thickness before
[11494]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))
[11403]215      !
[11395]216      ! get ktbl
[12340]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
[11395]226      !
[11541]227   END SUBROUTINE isf_tbl_kbot
[11395]228      !
[11541]229   SUBROUTINE isf_tbl_ktop(pdep, ktop)
[11403]230      !!--------------------------------------------------------------------
[11541]231      !!                  ***  ROUTINE isf_tbl_top  ***
[11395]232      !!
233      !! ** Purpose : compute top level of the isf top boundary layer in case of an ice shelf parametrisation
234      !!
235      !!-------------------------- OUT -------------------------------------
[11403]236      INTEGER,  DIMENSION(jpi,jpj), INTENT(  out) :: ktop        ! top level affected by the ice shelf parametrisation
[11395]237      !!-------------------------- IN  -------------------------------------
[11876]238      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pdep        ! top depth of the parametrisation influence
[11403]239      !!--------------------------------------------------------------------
[11395]240      INTEGER :: ji,jj
241      INTEGER :: ikt
[11403]242      !!--------------------------------------------------------------------
243      !
[11876]244      ! if we need to recompute the top level at every time stepcompute top level (z*, z~)
[12068]245      ! in case of weak ht variation we can assume the top level of htbl to be constant
[11876]246      ! => only done using gdepw_0
[11395]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
[11403]250      !
[12340]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
[11395]260      !
[11541]261   END SUBROUTINE isf_tbl_ktop
[11395]262
263END MODULE isftbl
Note: See TracBrowser for help on using the repository browser.