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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynspg.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: 11.9 KB
Line 
1MODULE dynspg
2   !!======================================================================
3   !!                       ***  MODULE  dynspg  ***
4   !! Ocean dynamics:  surface pressure gradient control
5   !!======================================================================
6   !! History :  1.0  ! 2005-12  (C. Talandier, G. Madec, V. Garnier)  Original code
7   !!            3.2  ! 2009-07  (R. Benshila)  Suppression of rigid-lid option
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   dyn_spg     : update the dynamics trend with surface pressure gradient
12   !!   dyn_spg_init: initialization, namelist read, and parameters control
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers variables
15   USE dom_oce        ! ocean space and time domain variables
16   USE c1d            ! 1D vertical configuration
17   USE phycst         ! physical constants
18   USE sbc_oce        ! surface boundary condition: ocean
19   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
20   USE sbcapr         ! surface boundary condition: atmospheric pressure
21   USE dynspg_exp     ! surface pressure gradient     (dyn_spg_exp routine)
22   USE dynspg_ts      ! surface pressure gradient     (dyn_spg_ts  routine)
23   USE tide_mod       !
24   USE trd_oce        ! trends: ocean variables
25   USE trddyn         ! trend manager: dynamics
26   !
27   USE prtctl         ! Print control                     (prt_ctl routine)
28   USE in_out_manager ! I/O manager
29   USE lib_mpp        ! MPP library
30   USE timing         ! Timing
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   dyn_spg        ! routine called by step module
36   PUBLIC   dyn_spg_init   ! routine called by opa module
37
38   INTEGER ::   nspg = 0   ! type of surface pressure gradient scheme defined from lk_dynspg_...
39
40   !                       ! Parameter to control the surface pressure gradient scheme
41   INTEGER, PARAMETER ::   np_TS  = 1   ! split-explicit time stepping (Time-Splitting)
42   INTEGER, PARAMETER ::   np_EXP = 0   !       explicit time stepping
43   INTEGER, PARAMETER ::   np_NO  =-1   ! no surface pressure gradient, no scheme
44   !
45   REAL(wp) ::   zt0step !   Time of day at the beginning of the time step
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49#  include "do_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
58      !!----------------------------------------------------------------------
59      !!                  ***  ROUTINE dyn_spg  ***
60      !!
61      !! ** Purpose :   compute surface pressure gradient including the
62      !!              atmospheric pressure forcing (ln_apr_dyn=T).
63      !!
64      !! ** Method  :   Two schemes:
65      !!              - explicit       : the spg is evaluated at now
66      !!              - split-explicit : a time splitting technique is used
67      !!
68      !!              ln_apr_dyn=T : the atmospheric pressure forcing is applied
69      !!             as the gradient of the inverse barometer ssh:
70      !!                apgu = - 1/rau0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
71      !!                apgv = - 1/rau0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
72      !!             Note that as all external forcing a time averaging over a two rdt
73      !!             period is used to prevent the divergence of odd and even time step.
74      !!----------------------------------------------------------------------
75      INTEGER                             , INTENT( in )  ::  kt                  ! ocean time-step index
76      INTEGER                             , INTENT( in )  ::  Kbb, Kmm, Krhs, Kaa ! ocean time level indices
77      REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) ::  puu, pvv            ! ocean velocities and RHS of momentum equation
78      REAL(wp), DIMENSION(jpi,jpj,jpt)    , INTENT(inout) ::  pssh, puu_b, pvv_b  ! SSH and barotropic velocities at main time levels
79      !
80      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
81      REAL(wp) ::   z2dt, zg_2, zintp, zgrau0r, zld   ! local scalars
82      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zpice
83      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv
84      !!----------------------------------------------------------------------
85      !
86      IF( ln_timing )   CALL timing_start('dyn_spg')
87      !
88      IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends
89         ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 
90         ztrdu(:,:,:) = puu(:,:,:,Krhs)
91         ztrdv(:,:,:) = pvv(:,:,:,Krhs)
92      ENDIF
93      !
94      IF(      ln_apr_dyn                                                &   ! atmos. pressure
95         .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) )   &   ! tide potential (no time slitting)
96         .OR.  ln_ice_embd ) THEN                                            ! embedded sea-ice
97         !
98         DO_2D_00_00
99            spgu(ji,jj) = 0._wp
100            spgv(ji,jj) = 0._wp
101         END_2D
102         !
103         IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN   !==  Atmospheric pressure gradient (added later in time-split case) ==!
104            zg_2 = grav * 0.5
105            DO_2D_00_00
106               spgu(ji,jj) = spgu(ji,jj) + zg_2 * (  ssh_ib (ji+1,jj) - ssh_ib (ji,jj)    &
107                  &                                + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
108               spgv(ji,jj) = spgv(ji,jj) + zg_2 * (  ssh_ib (ji,jj+1) - ssh_ib (ji,jj)    &
109                  &                                + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
110            END_2D
111         ENDIF
112         !
113         !                                    !==  tide potential forcing term  ==!
114         IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide )  ) THEN   ! N.B. added directly at sub-time-step in ts-case
115            !
116            ! Update tide potential at the beginning of current time step
117            zt0step = REAL(nsec_day, wp)-0.5_wp*rdt
118            CALL upd_tide(zt0step, Kmm)
119            !
120            DO_2D_00_00
121               spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
122               spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
123            END_2D
124            !
125            IF (ln_scal_load) THEN
126               zld = rn_scal_load * grav
127               DO_2D_00_00
128                  spgu(ji,jj) = spgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
129                  spgv(ji,jj) = spgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
130               END_2D
131            ENDIF
132         ENDIF
133         !
134         IF( ln_ice_embd ) THEN              !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
135            ALLOCATE( zpice(jpi,jpj) )
136            zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
137            zgrau0r     = - grav * r1_rau0
138            zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrau0r
139            DO_2D_00_00
140               spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
141               spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
142            END_2D
143            DEALLOCATE( zpice )         
144         ENDIF
145         !
146         DO_3D_00_00( 1, jpkm1 )
147            puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj)
148            pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + spgv(ji,jj)
149         END_3D
150         !
151!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
152         !   
153      ENDIF
154      !
155      SELECT CASE ( nspg )                   !== surface pressure gradient computed and add to the general trend ==!
156      CASE ( np_EXP )   ;   CALL dyn_spg_exp( kt,      Kmm,       puu, pvv, Krhs )                    ! explicit
157      CASE ( np_TS  )   ;   CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting
158      END SELECT
159      !                   
160      IF( l_trddyn )   THEN                  ! save the surface pressure gradient trends for further diagnostics
161         ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
162         ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
163         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt, Kmm )
164         DEALLOCATE( ztrdu , ztrdv ) 
165      ENDIF
166      !                                      ! print mean trends (used for debugging)
167      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' spg  - Ua: ', mask1=umask, &
168         &                                  tab3d_2=pvv(:,:,:,Krhs), clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
169      !
170      IF( ln_timing )   CALL timing_stop('dyn_spg')
171      !
172   END SUBROUTINE dyn_spg
173
174
175   SUBROUTINE dyn_spg_init
176      !!---------------------------------------------------------------------
177      !!                  ***  ROUTINE dyn_spg_init  ***
178      !!               
179      !! ** Purpose :   Control the consistency between namelist options for
180      !!              surface pressure gradient schemes
181      !!----------------------------------------------------------------------
182      INTEGER ::   ioptio, ios   ! local integers
183      !
184      NAMELIST/namdyn_spg/ ln_dynspg_exp       , ln_dynspg_ts,   &
185      &                    ln_bt_fw, ln_bt_av  , ln_bt_auto  ,   &
186      &                    nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
187      !!----------------------------------------------------------------------
188      !
189      IF(lwp) THEN
190         WRITE(numout,*)
191         WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme'
192         WRITE(numout,*) '~~~~~~~~~~~~'
193      ENDIF
194      !
195      READ  ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901)
196901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_spg in reference namelist' )
197      !
198      READ  ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 )
199902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist' )
200      IF(lwm) WRITE ( numond, namdyn_spg )
201      !
202      IF(lwp) THEN             ! Namelist print
203         WRITE(numout,*) '   Namelist : namdyn_spg                    '
204         WRITE(numout,*) '      Explicit free surface                  ln_dynspg_exp = ', ln_dynspg_exp
205         WRITE(numout,*) '      Free surface with time splitting       ln_dynspg_ts  = ', ln_dynspg_ts
206      ENDIF
207      !                          ! Control of surface pressure gradient scheme options
208                                     nspg =  np_NO    ;   ioptio = 0
209      IF( ln_dynspg_exp ) THEN   ;   nspg =  np_EXP   ;   ioptio = ioptio + 1   ;   ENDIF
210      IF( ln_dynspg_ts  ) THEN   ;   nspg =  np_TS    ;   ioptio = ioptio + 1   ;   ENDIF
211      !
212      IF( ioptio  > 1 )   CALL ctl_stop( 'Choose only one surface pressure gradient scheme' )
213      IF( ioptio == 0 )   CALL ctl_warn( 'NO surface pressure gradient trend in momentum Eqs.' )
214      IF( ln_dynspg_exp .AND. ln_isfcav )   &
215           &   CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' )
216      !
217      IF(lwp) THEN
218         WRITE(numout,*)
219         IF( nspg == np_EXP )   WRITE(numout,*) '   ==>>>   explicit free surface'
220         IF( nspg == np_TS  )   WRITE(numout,*) '   ==>>>   free surface with time splitting scheme'
221         IF( nspg == np_NO  )   WRITE(numout,*) '   ==>>>   No surface surface pressure gradient trend in momentum Eqs.'
222      ENDIF
223      !
224      IF( nspg == np_TS ) THEN   ! split-explicit scheme initialisation
225         CALL dyn_spg_ts_init          ! do it first: set nn_baro used to allocate some arrays later on
226      ENDIF
227      !
228   END SUBROUTINE dyn_spg_init
229
230  !!======================================================================
231END MODULE dynspg
Note: See TracBrowser for help on using the repository browser.