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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/USR/usrdef_hgr.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: 9.5 KB
Line 
1MODULE usrdef_hgr
2   !!======================================================================
3   !!                     ***  MODULE usrdef_hgr   ***
4   !!
5   !!                     ===  GYRE configuration  ===
6   !!
7   !! User defined :   mesh and Coriolis parameter of a user configuration
8   !!======================================================================
9   !! History :  4.0 ! 2016-03  (S. Flavoni)
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   usr_def_hgr   : initialize the horizontal mesh
14   !!----------------------------------------------------------------------
15   USE dom_oce  , ONLY: nimpp, njmpp       ! ocean space and time domain
16   USE par_oce        ! ocean space and time domain
17   USE phycst         ! physical constants
18   USE usrdef_nam     !
19   !
20   USE in_out_manager ! I/O manager
21   USE lib_mpp        ! MPP library
22   
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   usr_def_hgr   ! called in domhgr.F90
27
28   !! * Substitutions
29#  include "do_loop_substitute.h90"
30   !!----------------------------------------------------------------------
31   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
32   !! $Id$
33   !! Software governed by the CeCILL license (see ./LICENSE)
34   !!----------------------------------------------------------------------
35CONTAINS
36
37   SUBROUTINE usr_def_hgr( plamt , plamu , plamv  , plamf  ,   &   ! geographic position (required)
38      &                    pphit , pphiu , pphiv  , pphif  ,   &   !
39      &                    kff   , pff_f , pff_t  ,            &   ! Coriolis parameter  (if domain not on the sphere)
40      &                    pe1t  , pe1u  , pe1v   , pe1f   ,   &   ! scale factors       (required)
41      &                    pe2t  , pe2u  , pe2v   , pe2f   ,   &   !
42      &                    ke1e2u_v      , pe1e2u , pe1e2v     )   ! u- & v-surfaces (if gridsize reduction is used in strait(s))
43      !!----------------------------------------------------------------------
44      !!                  ***  ROUTINE usr_def_hgr  ***
45      !!
46      !! ** Purpose :   user defined mesh and Coriolis parameter
47      !!
48      !! ** Method  :   set all intent(out) argument to a proper value
49      !!
50      !!                Here GYRE configuration :
51      !!          Rectangular mid-latitude domain
52      !!          - with axes rotated by 45 degrees
53      !!          - a constant horizontal resolution of 106 km
54      !!          - on a beta-plane
55      !!
56      !! ** Action  : - define longitude & latitude of t-, u-, v- and f-points (in degrees)
57      !!              - define coriolis parameter at f-point if the domain in not on the sphere (on beta-plane)
58      !!              - define i- & j-scale factors at t-, u-, v- and f-points (in meters)
59      !!              - define u- & v-surfaces (if gridsize reduction is used in some straits) (in m2)
60      !!----------------------------------------------------------------------
61      REAL(wp), DIMENSION(:,:), INTENT(out) ::   plamt, plamu, plamv, plamf   ! longitude outputs                     [degrees]
62      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pphit, pphiu, pphiv, pphif   ! latitude outputs                      [degrees]
63      INTEGER                 , INTENT(out) ::   kff                          ! =1 Coriolis parameter computed here, =0 otherwise
64      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pff_f, pff_t                 ! Coriolis factor at f-point                [1/s]
65      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1t, pe1u, pe1v, pe1f       ! i-scale factors                             [m]
66      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe2t, pe2u, pe2v, pe2f       ! j-scale factors                             [m]
67      INTEGER                 , INTENT(out) ::   ke1e2u_v                     ! =1 u- & v-surfaces computed here, =0 otherwise
68      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pe1e2u, pe1e2v               ! u- & v-surfaces (if reduction in strait)   [m2]
69      !
70      INTEGER  ::   ji, jj               ! dummy loop indices
71      REAL(wp) ::   zlam1, zlam0, zcos_alpha, zim1 , zjm1 , ze1  , ze1deg, zf0 ! local scalars
72      REAL(wp) ::   zphi1, zphi0, zsin_alpha, zim05, zjm05, zbeta, znorme      !   -      -
73      !!-------------------------------------------------------------------------------
74      !
75      !     !==  beta-plane with regular grid-spacing and rotated domain ==!  (GYRE configuration)
76      !
77      IF(lwp) WRITE(numout,*)
78      IF(lwp) WRITE(numout,*) 'usr_def_hgr : GYRE configuration (beta-plane with rotated regular grid-spacing)'
79      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
80      !
81      !
82      !                       !==  grid point position  ==!
83      !
84      zlam1 = -85._wp                           ! position of gridpoint (i,j) = (1,jpjglo)
85      zphi1 =  29._wp
86      !
87      ze1 = 106000._wp / REAL( nn_GYRE , wp )   ! gridspacing in meters
88      !
89      zsin_alpha = - SQRT( 2._wp ) * 0.5_wp     ! angle: 45 degrees
90      zcos_alpha =   SQRT( 2._wp ) * 0.5_wp
91      ze1deg = ze1 / (ra * rad)
92      zlam0 = zlam1 + zcos_alpha * ze1deg * REAL( jpjglo-2 , wp )
93      zphi0 = zphi1 + zsin_alpha * ze1deg * REAL( jpjglo-2 , wp )
94
95#if defined key_agrif
96      ! ! Upper left longitude and latitude from parent:
97      IF (.NOT.Agrif_root()) THEN
98         zlam0 = zlam1 + Agrif_irhox() * REAL(Agrif_Parent(jpjglo)-2 , wp) * ze1deg * zcos_alpha  &
99                   &   + ( Agrif_Ix()*Agrif_irhox()-(0.5_wp+nbghostcells)) * ze1deg * zcos_alpha  &
100                   &   + ( Agrif_Iy()*Agrif_irhoy()-(0.5_wp+nbghostcells)) * ze1deg * zsin_alpha
101         zphi0 = zphi1 + Agrif_irhoy() * REAL(Agrif_Parent(jpjglo)-2 , wp) * ze1deg * zsin_alpha  &
102                   &   - ( Agrif_Ix()*Agrif_irhox()-nbghostcells )         * ze1deg * zsin_alpha  &
103                   &   + ( Agrif_Iy()*Agrif_irhoy()-nbghostcells )         * ze1deg * zcos_alpha
104      ENDIF 
105#endif
106      !   
107      IF( ln_bench ) THEN     ! benchmark: forced the resolution to be 106 km
108         ze1 = 106000._wp     ! but keep (lat,lon) at the right nn_GYRE resolution
109         CALL ctl_warn( ' GYRE used as Benchmark: e1=e2=106km, no need to adjust rdt, ahm,aht ' )
110      ENDIF
111      IF( nprint==1 .AND. lwp )   THEN
112         WRITE(numout,*) 'ze1', ze1, 'cosalpha', zcos_alpha, 'sinalpha', zsin_alpha
113         WRITE(numout,*) 'ze1deg', ze1deg, 'zlam0', zlam0, 'zphi0', zphi0
114      ENDIF
115      !   
116      DO_2D_11_11
117         zim1 = REAL( ji + nimpp - 1 ) - 1.   ;   zim05 = REAL( ji + nimpp - 1 ) - 1.5 
118         zjm1 = REAL( jj + njmpp - 1 ) - 1.   ;   zjm05 = REAL( jj + njmpp - 1 ) - 1.5 
119         !   
120         !glamt(i,j) longitude at T-point
121         !gphit(i,j) latitude at T-point 
122         plamt(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
123         pphit(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
124         !   
125         !glamu(i,j) longitude at U-point
126         !gphiu(i,j) latitude at U-point
127         plamu(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm05 * ze1deg * zsin_alpha
128         pphiu(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm05 * ze1deg * zcos_alpha
129         !   
130         !glamv(i,j) longitude at V-point
131         !gphiv(i,j) latitude at V-point
132         plamv(ji,jj) = zlam0 + zim05 * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
133         pphiv(ji,jj) = zphi0 - zim05 * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
134         !
135         !glamf(i,j) longitude at F-point
136         !gphif(i,j) latitude at F-point
137         plamf(ji,jj) = zlam0 + zim1  * ze1deg * zcos_alpha + zjm1  * ze1deg * zsin_alpha
138         pphif(ji,jj) = zphi0 - zim1  * ze1deg * zsin_alpha + zjm1  * ze1deg * zcos_alpha
139      END_2D
140      !
141      !                       !== Horizontal scale factors ==! (in meters)
142      !                     
143      !                                         ! constant grid spacing
144      pe1t(:,:) =  ze1     ;      pe2t(:,:) = ze1
145      pe1u(:,:) =  ze1     ;      pe2u(:,:) = ze1
146      pe1v(:,:) =  ze1     ;      pe2v(:,:) = ze1
147      pe1f(:,:) =  ze1     ;      pe2f(:,:) = ze1
148      !
149      !                                         ! NO reduction of grid size in some straits
150      ke1e2u_v = 0                              !    ==>> u_ & v_surfaces will be computed in dom_ghr routine
151      pe1e2u(:,:) = 0._wp                       !    CAUTION: set to zero to avoid error with some compilers that
152      pe1e2v(:,:) = 0._wp                       !             require an initialization of INTENT(out) arguments
153      !
154      !
155      !                       !==  Coriolis parameter  ==!
156      kff = 1                                            !  indicate not to compute ff afterward
157      !
158      zbeta = 2. * omega * COS( rad * zphi1 ) / ra       ! beta at latitude zphi1
159      !SF we overwrite zphi0 (south point in latitude) used just above to define pphif (value of zphi0=15.5190567531966)
160      !SF for computation of Coriolis we keep the parameter of Hazeleger, W., and S. S. Drijfhout, JPO 1998.
161      zphi0 = 15._wp                                     !  latitude of the most southern grid point 
162      zf0   = 2. * omega * SIN( rad * zphi0 )            !  compute f0 1st point south
163      !
164      pff_f(:,:) = ( zf0 + zbeta * ABS( pphif(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
165      pff_t(:,:) = ( zf0 + zbeta * ABS( pphit(:,:) - zphi0 ) * rad * ra ) ! f = f0 +beta* y ( y=0 at south)
166      !
167      IF(lwp) WRITE(numout,*) '                           beta-plane used. beta = ', zbeta, ' 1/(s.m)'
168      !
169   END SUBROUTINE usr_def_hgr
170
171   !!======================================================================
172END MODULE usrdef_hgr
Note: See TracBrowser for help on using the repository browser.