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

source: NEMO/branches/UKMO/BPC_miniapp/Master/dynhpgs.F90 @ 10888

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

Ticket #2197 - added the Master directory with code and makefile
File size: 16.1 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!$OMP PARALLEL DO PRIVATE(zcoef1)
84      DO jj = 2, jpjm1
85         DO ji = fs_2, fs_jpim1   ! vector opt.
86            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
87            ! hydrostatic pressure gradient
88            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
89            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
90            ! add to the general momentum trend
91            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
92            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
93         END DO
94      END DO
95
96      !
97      ! interior value (2=<jk=<jpkm1)
98!$OMP PARALLEL DO PRIVATE(zcoef1)
99      DO jk = 2, jpkm1
100         DO jj = 2, jpjm1
101            DO ji = fs_2, fs_jpim1   ! vector opt.
102               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
103               ! hydrostatic pressure gradient
104               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
105                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    &
106                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
107
108               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
109                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    &
110                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
111               ! add to the general momentum trend
112               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
113               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
114            END DO
115         END DO
116      END DO
117      !
118!$ACC END KERNELS
119      !hpg_zco_time = hpg_zco_time + (TIMER() - et) ! Moved time up the call tree
120
121   END SUBROUTINE hpg_zco
122
123
124   SUBROUTINE hpg_zps( grav, r1_e1u, r1_e2v, mbku, mbkv, gru, grv, e3w_n, rhd, ua, va, zhpi, zhpj )
125      !!---------------------------------------------------------------------
126      !!                 ***  ROUTINE hpg_zps  ***
127      !!
128      !! ** Method  :   z-coordinate plus partial steps case.  blahblah...
129      !!
130      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
131      !!----------------------------------------------------------------------
132
133      REAL(wp), INTENT(in)::   grav
134      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj) ::      r1_e1u, r1_e2v, gru, grv
135      INTEGER, INTENT(in),     DIMENSION(jpi,jpj) ::      mbku, mbkv
136      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj,jpk) ::  e3w_n, rhd
137      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  ua, va 
138      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
139
140      !!----------------------------------------------------------------------
141      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
142      INTEGER  ::   iku, ikv                         ! temporary integers
143      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
144      REAL(wp) :: et
145
146      !!----------------------------------------------------------------------
147      !
148! suppressed by MJB 27/10/2018 - should be moved to a control routine
149!      IF( kt == nit000 ) THEN
150!         IF(lwp) WRITE(numout,*)
151!         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
152!         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
153!      ENDIF
154
155      ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level
156!jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv )
157      et = TIMER()
158!$ACC KERNELS
159
160      ! Local constant initialization
161      zcoef0 = - grav * 0.5_wp
162
163
164      !  Surface value (also valid in partial step case)
165!$OMP PARALLEL DO PRIVATE(zcoef1)
166      DO jj = 2, jpjm1
167         DO ji = fs_2, fs_jpim1   ! vector opt.
168            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
169            ! hydrostatic pressure gradient
170            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
171            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
172            ! add to the general momentum trend
173            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
174            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
175         END DO
176      END DO
177
178      ! interior value (2=<jk=<jpkm1)
179!$OMP PARALLEL DO PRIVATE(zcoef1)
180      DO jk = 2, jpkm1
181         DO jj = 2, jpjm1
182            DO ji = fs_2, fs_jpim1   ! vector opt.
183               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
184               ! hydrostatic pressure gradient
185               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
186                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
187                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
188
189               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
190                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
191                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
192               ! add to the general momentum trend
193               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
194               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
195            END DO
196         END DO
197      END DO
198
199      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
200!$OMP PARALLEL DO PRIVATE(zcoef1,zcoef2,zcoef3,iku,ikv)
201      DO jj = 2, jpjm1
202         DO ji = 2, jpim1
203            iku = mbku(ji,jj)
204            ikv = mbkv(ji,jj)
205            zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) )
206            zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) )
207            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
208               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
209               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
210                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
211               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
212            ENDIF
213            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
214               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
215               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
216                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
217               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
218            ENDIF
219         END DO
220      END DO
221!$ACC END KERNELS
222      !
223      !hpg_zps_time = hpg_zps_time + (TIMER() - et) ! Moved timer up call tree
224
225   END SUBROUTINE hpg_zps
226
227
228   SUBROUTINE hpg_sco( grav, r1_e1u, r1_e2v, e3w_n, rhd, gde3w_n, ln_linssh, ua, va, zhpi, zhpj)
229      !!---------------------------------------------------------------------
230      !!                  ***  ROUTINE hpg_sco  ***
231      !!
232      !! ** Method  :   s-coordinate case. Jacobian scheme.
233      !!      The now hydrostatic pressure gradient at a given level, jk,
234      !!      is computed by taking the vertical integral of the in-situ
235      !!      density gradient along the model level from the suface to that
236      !!      level. s-coordinates (ln_sco): a corrective term is added
237      !!      to the horizontal pressure gradient :
238      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
239      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
240      !!      add it to the general momentum trend (ua,va).
241      !!         ua = ua - 1/e1u * zhpi
242      !!         va = va - 1/e2v * zhpj
243      !!
244      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
245      !!----------------------------------------------------------------------
246      REAL(wp), INTENT(in)::   grav
247      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj) ::      r1_e1u, r1_e2v
248      REAL(wp), INTENT(in),    DIMENSION(jpi,jpj,jpk) ::  e3w_n, rhd, gde3w_n
249      LOGICAL,  INTENT(in)::   ln_linssh
250      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::  ua, va 
251      REAL(wp), INTENT(INOUT), DIMENSION(jpi,jpj,jpk) ::   zhpi, zhpj
252
253      !!----------------------------------------------------------------------
254      !!
255      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
256      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
257      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
258      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zcpx, zcpy   !W/D pressure filter
259      REAL(wp) :: et
260      !!----------------------------------------------------------------------
261      !
262! suppressed by MJB 27/10/2018 - should be moved to a control routine
263!      IF( kt == nit000 ) THEN
264!         IF(lwp) WRITE(numout,*)
265!         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
266!         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
267!      ENDIF
268      !
269      et = TIMER()
270      zcoef0 = - grav * 0.5_wp
271      IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly
272      ELSE                    ;   znad = 1._wp         ! Variable volume: density
273      ENDIF
274      !
275!$ACC KERNELS
276      ! Surface value
277!$OMP PARALLEL
278!$OMP DO PRIVATE(zuap,zvap)
279      DO jj = 2, jpjm1
280         DO ji = fs_2, fs_jpim1   ! vector opt.
281            ! hydrostatic pressure gradient along s-surfaces
282            zhpi(ji,jj,1) = zcoef0 * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )    &
283               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
284            zhpj(ji,jj,1) = zcoef0 * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )    &
285               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
286            ! s-coordinate pressure gradient correction
287            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
288               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
289            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
290               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
291            !
292            ! add to the general momentum trend
293            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
294            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
295         END DO
296      END DO
297
298      ! interior value (2=<jk=<jpkm1)
299!$OMP DO PRIVATE(zuap,zvap)
300      DO jk = 2, jpkm1
301         DO jj = 2, jpjm1
302            DO ji = fs_2, fs_jpim1   ! vector opt.
303               ! hydrostatic pressure gradient along s-surfaces
304               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
305                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
306                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
307               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
308                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
309                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
310               ! s-coordinate pressure gradient correction
311               zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
312                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)
313               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
314                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
315               !
316               ! add to the general momentum trend
317               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
318               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
319            END DO
320         END DO
321      END DO
322!$OMP END PARALLEL
323!$ACC END KERNELS
324      !
325      !hpg_sco_time = hpg_sco_time + (TIMER() - et) ! Moved timer up the call tree
326   END SUBROUTINE hpg_sco
327
328   !!======================================================================
329END MODULE dynhpgs
330
Note: See TracBrowser for help on using the repository browser.