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.
dynhpgs.F90 in NEMO/branches/UKMO/BPC_miniapp/OpenACC_managed – NEMO

source: NEMO/branches/UKMO/BPC_miniapp/OpenACC_managed/dynhpgs.F90 @ 10838

Last change on this file since 10838 was 10838, checked in by wayne_gaudin, 5 years ago

Ticket #2197 - extracted versions added
File size: 15.9 KB
Line 
1MODULE dynhpgs
2   !!======================================================================
3   !!                       ***  MODULE  dynhpgs  ***
4   !! Ocean dynamics:  hydrostatic pressure gradient trend ; science routines
5   !!======================================================================
6   !! History :  OPA  !  1987-09  (P. Andrich, M.-A. Foujols)  hpg_zco: Original code
7   !!            5.0  !  1991-11  (G. Madec)
8   !!            7.0  !  1996-01  (G. Madec)  hpg_sco: Original code for s-coordinates
9   !!            8.0  !  1997-05  (G. Madec)  split dynber into dynkeg and dynhpg
10   !!            8.5  !  2002-07  (G. Madec)  F90: Free form and module
11   !!            8.5  !  2002-08  (A. Bozec)  hpg_zps: Original code
12   !!   NEMO     1.0  !  2005-10  (A. Beckmann, B.W. An)  various s-coordinate options
13   !!                 !         Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot
14   !!             -   !  2005-11  (G. Madec) style & small optimisation
15   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
16   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates
17   !!                 !           (A. Coward) suppression of hel, wdj and rot options
18   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!       hpg_zco  : z-coordinate scheme
23   !!       hpg_zps  : z-coordinate plus partial steps (interpolation)
24   !!       hpg_sco  : s-coordinate (standard jacobian formulation)
25   !!----------------------------------------------------------------------
26   !
27USE par_kind
28USE len_oce
29USE in_out_manager
30
31   IMPLICIT NONE
32   !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
36   !! $Id: dynhpg.F90 9490 2018-04-23 08:44:07Z gm $
37   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41
42   SUBROUTINE hpg_zco( grav, r1_e1u, r1_e2v, e3w_n, rhd, ua, va, zhpi, zhpj)
43      !!---------------------------------------------------------------------
44      !!                  ***  ROUTINE hpg_zco  ***
45      !!
46      !! ** Method  :   z-coordinate case, levels are horizontal surfaces.
47      !!      The now hydrostatic pressure gradient at a given level, jk,
48      !!      is computed by taking the vertical integral of the in-situ
49      !!      density gradient along the model level from the suface to that
50      !!      level:    zhpi = grav .....
51      !!                zhpj = grav .....
52      !!      add it to the general momentum trend (ua,va).
53      !!            ua = ua - 1/e1u * zhpi
54      !!            va = va - 1/e2v * zhpj
55      !!
56      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
57      !!----------------------------------------------------------------------
58      REAL(wp), INTENT(in) ::                             grav 
59      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj) ::      r1_e1u, r1_e2v
60      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj,jpk) ::  e3w_n, rhd
61      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  ua, va 
62      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
63      REAL(wp) :: et
64
65      !!----------------------------------------------------------------------
66      !
67      INTEGER  ::   ji, jj, jk       ! dummy loop indices
68      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
69      !!----------------------------------------------------------------------
70      !
71! suppressed by MJB 27/10/2018 - should be moved to a control routine
72!      IF( kt == nit000 ) THEN
73!         IF(lwp) WRITE(numout,*)
74!         IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
75!         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
76!      ENDIF
77      et = TIMER()
78
79!$ACC KERNELS
80      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
81
82      ! Surface value
83      DO jj = 2, jpjm1
84         DO ji = fs_2, fs_jpim1   ! vector opt.
85            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
86            ! hydrostatic pressure gradient
87            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
88            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
89            ! add to the general momentum trend
90            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
91            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
92         END DO
93      END DO
94
95      !
96      ! interior value (2=<jk=<jpkm1)
97      DO jk = 2, jpkm1
98         DO jj = 2, jpjm1
99            DO ji = fs_2, fs_jpim1   ! vector opt.
100               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
101               ! hydrostatic pressure gradient
102               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
103                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    &
104                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
105
106               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
107                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    &
108                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
109               ! add to the general momentum trend
110               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
111               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
112            END DO
113         END DO
114      END DO
115      !
116!$ACC END KERNELS
117      !hpg_zco_time = hpg_zco_time + (TIMER() - et) ! Moved time up the call tree
118
119   END SUBROUTINE hpg_zco
120
121
122   SUBROUTINE hpg_zps( grav, r1_e1u, r1_e2v, mbku, mbkv, gru, grv, e3w_n, rhd, ua, va, zhpi, zhpj )
123      !!---------------------------------------------------------------------
124      !!                 ***  ROUTINE hpg_zps  ***
125      !!
126      !! ** Method  :   z-coordinate plus partial steps case.  blahblah...
127      !!
128      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
129      !!----------------------------------------------------------------------
130
131      REAL(wp), INTENT(in)::   grav
132      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj) ::      r1_e1u, r1_e2v, gru, grv
133      INTEGER, INTENT(in),     DIMENSION(jpi,jpj) ::      mbku, mbkv
134      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj,jpk) ::  e3w_n, rhd
135      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  ua, va 
136      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
137
138      !!----------------------------------------------------------------------
139      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
140      INTEGER  ::   iku, ikv                         ! temporary integers
141      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
142      REAL(wp) :: et
143
144      !!----------------------------------------------------------------------
145      !
146! suppressed by MJB 27/10/2018 - should be moved to a control routine
147!      IF( kt == nit000 ) THEN
148!         IF(lwp) WRITE(numout,*)
149!         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
150!         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
151!      ENDIF
152
153      ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level
154!jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv )
155      et = TIMER()
156!$ACC KERNELS
157
158      ! Local constant initialization
159      zcoef0 = - grav * 0.5_wp
160
161
162      !  Surface value (also valid in partial step case)
163      DO jj = 2, jpjm1
164         DO ji = fs_2, fs_jpim1   ! vector opt.
165            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
166            ! hydrostatic pressure gradient
167            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
168            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
169            ! add to the general momentum trend
170            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
171            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
172         END DO
173      END DO
174
175      ! interior value (2=<jk=<jpkm1)
176      DO jk = 2, jpkm1
177         DO jj = 2, jpjm1
178            DO ji = fs_2, fs_jpim1   ! vector opt.
179               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
180               ! hydrostatic pressure gradient
181               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
182                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
183                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
184
185               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
186                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
187                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
188               ! add to the general momentum trend
189               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
190               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
191            END DO
192         END DO
193      END DO
194
195      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
196      DO jj = 2, jpjm1
197         DO ji = 2, jpim1
198            iku = mbku(ji,jj)
199            ikv = mbkv(ji,jj)
200            zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) )
201            zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) )
202            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
203               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
204               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
205                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
206               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
207            ENDIF
208            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
209               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
210               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
211                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
212               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
213            ENDIF
214         END DO
215      END DO
216!$ACC END KERNELS
217      !
218      !hpg_zps_time = hpg_zps_time + (TIMER() - et) ! Moved timer up call tree
219
220   END SUBROUTINE hpg_zps
221
222
223   SUBROUTINE hpg_sco( grav, r1_e1u, r1_e2v, e3w_n, rhd, gde3w_n, ln_linssh, ua, va, zhpi, zhpj)
224      !!---------------------------------------------------------------------
225      !!                  ***  ROUTINE hpg_sco  ***
226      !!
227      !! ** Method  :   s-coordinate case. Jacobian scheme.
228      !!      The now hydrostatic pressure gradient at a given level, jk,
229      !!      is computed by taking the vertical integral of the in-situ
230      !!      density gradient along the model level from the suface to that
231      !!      level. s-coordinates (ln_sco): a corrective term is added
232      !!      to the horizontal pressure gradient :
233      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
234      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
235      !!      add it to the general momentum trend (ua,va).
236      !!         ua = ua - 1/e1u * zhpi
237      !!         va = va - 1/e2v * zhpj
238      !!
239      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
240      !!----------------------------------------------------------------------
241      REAL(wp), INTENT(in)::   grav
242      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj) ::      r1_e1u, r1_e2v
243      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj,jpk) ::  e3w_n, rhd, gde3w_n
244      LOGICAL,  INTENT(in)::   ln_linssh
245      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  ua, va 
246      REAL(wp), INTENT(INOUT), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
247
248      !!----------------------------------------------------------------------
249      !!
250      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
251      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
252      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
253      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
254      REAL(wp) :: et
255      !!----------------------------------------------------------------------
256      !
257! suppressed by MJB 27/10/2018 - should be moved to a control routine
258!      IF( kt == nit000 ) THEN
259!         IF(lwp) WRITE(numout,*)
260!         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
261!         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
262!      ENDIF
263      !
264      et = TIMER()
265      zcoef0 = - grav * 0.5_wp
266      IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly
267      ELSE                    ;   znad = 1._wp         ! Variable volume: density
268      ENDIF
269      !
270!$ACC KERNELS
271      ! Surface value
272      DO jj = 2, jpjm1
273         DO ji = fs_2, fs_jpim1   ! vector opt.
274            ! hydrostatic pressure gradient along s-surfaces
275            zhpi(ji,jj,1) = zcoef0 * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )    &
276               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
277            zhpj(ji,jj,1) = zcoef0 * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )    &
278               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
279            ! s-coordinate pressure gradient correction
280            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
281               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
282            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
283               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
284            !
285            ! add to the general momentum trend
286            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
287            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
288         END DO
289      END DO
290
291      ! interior value (2=<jk=<jpkm1)
292      DO jk = 2, jpkm1
293         DO jj = 2, jpjm1
294            DO ji = fs_2, fs_jpim1   ! vector opt.
295               ! hydrostatic pressure gradient along s-surfaces
296               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
297                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
298                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
299               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
300                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
301                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
302               ! s-coordinate pressure gradient correction
303               zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
304                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)
305               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
306                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
307               !
308               ! add to the general momentum trend
309               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
310               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
311            END DO
312         END DO
313      END DO
314!$ACC END KERNELS
315      !
316      !hpg_sco_time = hpg_sco_time + (TIMER() - et) ! Moved timer up the call tree
317   END SUBROUTINE hpg_sco
318
319   !!======================================================================
320END MODULE dynhpgs
321
Note: See TracBrowser for help on using the repository browser.