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/OpenMP – NEMO

source: NEMO/branches/UKMO/BPC_miniapp/OpenMP/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: 16.0 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      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
80
81      ! Surface value
82!$OMP PARALLEL DO PRIVATE(zcoef1)
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!$OMP PARALLEL DO PRIVATE(zcoef1)
98      DO jk = 2, jpkm1
99         DO jj = 2, jpjm1
100            DO ji = fs_2, fs_jpim1   ! vector opt.
101               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
102               ! hydrostatic pressure gradient
103               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
104                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    &
105                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
106
107               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
108                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    &
109                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
110               ! add to the general momentum trend
111               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
112               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
113            END DO
114         END DO
115      END DO
116      !
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
157      ! Local constant initialization
158      zcoef0 = - grav * 0.5_wp
159
160
161      !  Surface value (also valid in partial step case)
162!$OMP PARALLEL DO PRIVATE(zcoef1)
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!$OMP PARALLEL DO PRIVATE(zcoef1)
177      DO jk = 2, jpkm1
178         DO jj = 2, jpjm1
179            DO ji = fs_2, fs_jpim1   ! vector opt.
180               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
181               ! hydrostatic pressure gradient
182               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
183                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
184                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
185
186               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
187                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
188                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
189               ! add to the general momentum trend
190               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
191               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
192            END DO
193         END DO
194      END DO
195
196      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
197!$OMP PARALLEL DO PRIVATE(zcoef1,zcoef2,zcoef3,iku,ikv)
198      DO jj = 2, jpjm1
199         DO ji = 2, jpim1
200            iku = mbku(ji,jj)
201            ikv = mbkv(ji,jj)
202            zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) )
203            zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) )
204            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
205               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
206               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
207                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
208               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
209            ENDIF
210            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
211               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
212               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
213                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
214               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
215            ENDIF
216         END DO
217      END DO
218      !
219      !hpg_zps_time = hpg_zps_time + (TIMER() - et) ! Moved timer up call tree
220
221   END SUBROUTINE hpg_zps
222
223
224   SUBROUTINE hpg_sco( grav, r1_e1u, r1_e2v, e3w_n, rhd, gde3w_n, ln_linssh, ua, va, zhpi, zhpj)
225      !!---------------------------------------------------------------------
226      !!                  ***  ROUTINE hpg_sco  ***
227      !!
228      !! ** Method  :   s-coordinate case. Jacobian scheme.
229      !!      The now hydrostatic pressure gradient at a given level, jk,
230      !!      is computed by taking the vertical integral of the in-situ
231      !!      density gradient along the model level from the suface to that
232      !!      level. s-coordinates (ln_sco): a corrective term is added
233      !!      to the horizontal pressure gradient :
234      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
235      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
236      !!      add it to the general momentum trend (ua,va).
237      !!         ua = ua - 1/e1u * zhpi
238      !!         va = va - 1/e2v * zhpj
239      !!
240      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
241      !!----------------------------------------------------------------------
242      REAL(wp), INTENT(in)::   grav
243      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj) ::      r1_e1u, r1_e2v
244      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj,jpk) ::  e3w_n, rhd, gde3w_n
245      LOGICAL,  INTENT(in)::   ln_linssh
246      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  ua, va 
247      REAL(wp), INTENT(INOUT), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
248
249      !!----------------------------------------------------------------------
250      !!
251      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
252      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
253      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
254      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
255      REAL(wp) :: et
256      !!----------------------------------------------------------------------
257      !
258! suppressed by MJB 27/10/2018 - should be moved to a control routine
259!      IF( kt == nit000 ) THEN
260!         IF(lwp) WRITE(numout,*)
261!         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
262!         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
263!      ENDIF
264      !
265      et = TIMER()
266      zcoef0 = - grav * 0.5_wp
267      IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly
268      ELSE                    ;   znad = 1._wp         ! Variable volume: density
269      ENDIF
270      !
271      ! Surface value
272!$OMP PARALLEL
273!$OMP DO PRIVATE(zuap,zvap)
274      DO jj = 2, jpjm1
275         DO ji = fs_2, fs_jpim1   ! vector opt.
276            ! hydrostatic pressure gradient along s-surfaces
277            zhpi(ji,jj,1) = zcoef0 * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )    &
278               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
279            zhpj(ji,jj,1) = zcoef0 * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )    &
280               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
281            ! s-coordinate pressure gradient correction
282            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
283               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
284            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
285               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
286            !
287            ! add to the general momentum trend
288            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
289            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
290         END DO
291      END DO
292
293      ! interior value (2=<jk=<jpkm1)
294!$OMP DO PRIVATE(zuap,zvap)
295      DO jk = 2, jpkm1
296         DO jj = 2, jpjm1
297            DO ji = fs_2, fs_jpim1   ! vector opt.
298               ! hydrostatic pressure gradient along s-surfaces
299               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
300                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
301                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
302               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
303                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
304                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
305               ! s-coordinate pressure gradient correction
306               zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
307                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)
308               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
309                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
310               !
311               ! add to the general momentum trend
312               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
313               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
314            END DO
315         END DO
316      END DO
317!$OMP END PARALLEL
318      !
319      !hpg_sco_time = hpg_sco_time + (TIMER() - et) ! Moved timer up the call tree
320   END SUBROUTINE hpg_sco
321
322   !!======================================================================
323END MODULE dynhpgs
324
Note: See TracBrowser for help on using the repository browser.