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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/LDF/ldfc1d_c2d.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.
  • Property svn:keywords set to Id
File size: 7.9 KB
Line 
1MODULE ldfc1d_c2d
2   !!======================================================================
3   !!                    ***  MODULE  ldfc1d_c2d  ***
4   !! Ocean physics:  profile and horizontal shape of lateral eddy coefficients
5   !!=====================================================================
6   !! History :  3.7  ! 2013-12  (G. Madec)  restructuration/simplification of aht/aeiv specification,
7   !!                 !                      add velocity dependent coefficient and optional read in file
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   ldf_c1d       : ah reduced by 1/4 on the vertical (tanh profile, inflection at 300m)
12   !!   ldf_c2d       : ah = F(e1,e2) (laplacian or = F(e1^3,e2^3) (bilaplacian)
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers
15   USE dom_oce        ! ocean space and time domain
16   USE phycst         ! physical constants
17   !
18   USE in_out_manager ! I/O manager
19   USE lib_mpp        ! distribued memory computing library
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC   ldf_c1d   ! called by ldftra and ldfdyn modules
26   PUBLIC   ldf_c2d   ! called by ldftra and ldfdyn modules
27
28   REAL(wp) ::   r1_2  = 0.5_wp           ! =1/2
29   REAL(wp) ::   r1_4  = 0.25_wp          ! =1/4
30   REAL(wp) ::   r1_12 = 1._wp / 12._wp   ! =1/12
31 
32   !! * Substitutions
33#  include "do_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id$
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE ldf_c1d( cd_type, pahs1, pahs2, pah1, pah2 )
42      !!----------------------------------------------------------------------
43      !!                  ***  ROUTINE ldf_c1d  ***
44      !!             
45      !! ** Purpose :   1D eddy diffusivity/viscosity coefficients
46      !!
47      !! ** Method  :   1D eddy diffusivity coefficients F( depth )
48      !!                Reduction by zratio from surface to bottom
49      !!                hyperbolic tangent profile with inflection point
50      !!                at zh=500m and a width of zw=200m
51      !!
52      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
53      !!             DYN      pah1, pah2 defined at T- and F-points
54      !!----------------------------------------------------------------------
55      CHARACTER(len=3)                , INTENT(in   ) ::   cd_type        ! DYNamique or TRAcers
56      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in   ) ::   pahs1, pahs2   ! surface value of eddy coefficient   [m2/s]
57      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pah1 , pah2    ! eddy coefficient                    [m2/s]
58      !
59      INTEGER  ::   ji, jj, jk      ! dummy loop indices
60      REAL(wp) ::   zh, zc, zdep1   ! local scalars
61      REAL(wp) ::   zw    , zdep2   !   -      -
62      REAL(wp) ::   zratio          !   -      -
63      !!----------------------------------------------------------------------
64      !
65      IF(lwp) WRITE(numout,*)
66      IF(lwp) WRITE(numout,*) '   ldf_c1d : set a given profile to eddy mixing coefficients'
67      !
68      ! initialization of the profile
69      zratio = 0.25_wp           ! surface/bottom ratio
70      zh =  500._wp              ! depth    of the inflection point [m]
71      zw =  1._wp / 200._wp      ! width^-1     -        -      -   [1/m]
72      !                          ! associated coefficient           [-]
73      zc = ( 1._wp - zratio ) / ( 1._wp + TANH( zh * zw) )
74      !
75      !
76      SELECT CASE( cd_type )        ! point of calculation
77      !
78      CASE( 'DYN' )                     ! T- and F-points
79         DO jk = jpkm1, 1, -1                ! pah1 at T-point
80            pah1(:,:,jk) = pahs1(:,:) * (  zratio + zc * ( 1._wp + TANH( - ( gdept_0(:,:,jk) - zh ) * zw) )  )
81         END DO
82         DO_3DS_10_10( jpkm1, 1, -1 )
83            zdep2 = (  gdept_0(ji,jj+1,jk) + gdept_0(ji+1,jj+1,jk)   &
84               &     + gdept_0(ji,jj  ,jk) + gdept_0(ji+1,jj  ,jk)  ) * r1_4
85            pah2(ji,jj,jk) = pahs2(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  )
86         END_3D
87         CALL lbc_lnk( 'ldfc1d_c2d', pah2, 'F', 1. )   ! Lateral boundary conditions
88         !
89      CASE( 'TRA' )                     ! U- and V-points (zdep1 & 2 are an approximation in zps-coord.)
90         DO_3DS_10_10( jpkm1, 1, -1 )
91            zdep1 = (  gdept_0(ji,jj,jk) + gdept_0(ji+1,jj,jk)  ) * 0.5_wp
92            zdep2 = (  gdept_0(ji,jj,jk) + gdept_0(ji,jj+1,jk)  ) * 0.5_wp
93            pah1(ji,jj,jk) = pahs1(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep1 - zh ) * zw) )  )
94            pah2(ji,jj,jk) = pahs2(ji,jj) * (  zratio + zc * ( 1._wp + TANH( - ( zdep2 - zh ) * zw) )  )
95         END_3D
96         ! Lateral boundary conditions
97         CALL lbc_lnk_multi( 'ldfc1d_c2d', pah1, 'U', 1. , pah2, 'V', 1. )   
98         !
99      CASE DEFAULT                        ! error
100         CALL ctl_stop( 'ldf_c1d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' )
101      END SELECT
102      !
103   END SUBROUTINE ldf_c1d
104
105
106   SUBROUTINE ldf_c2d( cd_type, pUfac, knn, pah1, pah2 )
107      !!----------------------------------------------------------------------
108      !!                  ***  ROUTINE ldf_c2d  ***
109      !!             
110      !! ** Purpose :   2D eddy diffusivity/viscosity coefficients
111      !!
112      !! ** Method  :   2D eddy diffusivity coefficients F( e1 , e2 )
113      !!       laplacian   operator :   ah proportional to the scale factor      [m2/s]
114      !!       bilaplacian operator :   ah proportional to the (scale factor)^3  [m4/s]
115      !!       In both cases, pah0 is the maximum value reached by the coefficient
116      !!       at the Equator in case of e1=ra*rad= ~111km, not over the whole domain.
117      !!
118      !!   cd_type = TRA      pah1, pah2 defined at U- and V-points
119      !!             DYN      pah1, pah2 defined at T- and F-points
120      !!----------------------------------------------------------------------
121      CHARACTER(len=3)          , INTENT(in   ) ::   cd_type      ! DYNamique or TRAcers
122      REAL(wp)                  , INTENT(in   ) ::   pUfac        ! =1/2*Uc LAPlacian BiLaPlacian
123      INTEGER                   , INTENT(in   ) ::   knn          ! characteristic velocity   [m/s]
124      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   pah1, pah2   ! eddy coefficients         [m2/s or m4/s]
125      !
126      INTEGER  ::   ji, jj, jk   ! dummy loop indices
127      INTEGER  ::   inn          ! local integer
128      !!----------------------------------------------------------------------
129      !
130      IF(lwp) WRITE(numout,*)
131      IF(lwp) WRITE(numout,*) '   ldf_c2d :   aht = Ufac * max(e1,e2)   with Ufac = ', pUfac, ' m/s'
132      !
133      !
134      SELECT CASE( cd_type )        !==  surface values  ==!  (chosen grid point function of DYN or TRA)
135      !
136      CASE( 'DYN' )                       ! T- and F-points
137         DO_2D_11_11
138            pah1(ji,jj,1) = pUfac * MAX( e1t(ji,jj) , e2t(ji,jj) )**knn
139            pah2(ji,jj,1) = pUfac * MAX( e1f(ji,jj) , e2f(ji,jj) )**knn
140         END_2D
141      CASE( 'TRA' )                       ! U- and V-points
142         DO_2D_11_11
143            pah1(ji,jj,1) = pUfac * MAX( e1u(ji,jj), e2u(ji,jj) )**knn
144            pah2(ji,jj,1) = pUfac * MAX( e1v(ji,jj), e2v(ji,jj) )**knn
145         END_2D
146      CASE DEFAULT                        ! error
147         CALL ctl_stop( 'ldf_c2d: ', cd_type, ' Unknown, i.e. /= DYN or TRA' )
148      END SELECT
149      !                             !==  deeper values = surface one  ==!  (except jpk)
150      DO jk = 2, jpkm1
151         pah1(:,:,jk) = pah1(:,:,1)
152         pah2(:,:,jk) = pah2(:,:,1)
153      END DO
154      !
155   END SUBROUTINE ldf_c2d
156
157   !!======================================================================
158END MODULE ldfc1d_c2d
Note: See TracBrowser for help on using the repository browser.