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.
dynhpg.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago

Initial revision
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.0 KB
Line 
1MODULE dynhpg
2   !!======================================================================
3   !!                       ***  MODULE  dynhpg  ***
4   !! Ocean dynamics:  hydrostatic pressure gradient trend
5   !!======================================================================
6
7   !!----------------------------------------------------------------------
8   !!   dyn_hpg      : update the momentum trend with the horizontal
9   !!                  gradient of the hydrostatic pressure
10   !!
11   !!   default case : use of 3D work arrays (vector opt. available)
12   !!   key_s_coord       : s-coordinate
13   !!   key_partial_steps : z-coordinate with partial steps
14   !!   default key       : z-coordinate
15   !!----------------------------------------------------------------------
16   !! * Modules used
17   USE oce             ! ocean dynamics and tracers
18   USE dom_oce         ! ocean space and time domain
19   USE phycst          ! physical constants
20   USE in_out_manager  ! I/O manager
21   USE trddyn_oce     ! dynamics trends diagnostics variables
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Accessibility
27   PUBLIC dyn_hpg                ! routine called by step.F90
28
29#if defined key_autotasking
30   !!----------------------------------------------------------------------
31   !!   'key_autotasking' :                             j-k-i loop (j-slab)
32   !!----------------------------------------------------------------------
33   LOGICAL, PUBLIC ::   l_dyn_hpg_tsk = .TRUE.    ! ???
34   LOGICAL, PUBLIC ::   l_dyn_hpg     = .FALSE.   ! ???
35#else
36   !!----------------------------------------------------------------------
37   !!   default case :                             k-j-i loop (vector opt.)
38   !!----------------------------------------------------------------------   
39   LOGICAL, PUBLIC ::   l_dyn_hpg_tsk = .FALSE.   ! ???
40   LOGICAL, PUBLIC ::   l_dyn_hpg     = .TRUE.    ! ???
41#endif
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "vectopt_loop_substitute.h90"
46   !!----------------------------------------------------------------------
47   !!   OPA 9.0 , LODYC-IPSL  (2003)
48   !!----------------------------------------------------------------------
49
50CONTAINS
51
52#if defined key_s_coord
53   !!----------------------------------------------------------------------
54   !!   'key_s_coord' :                                        s-coordinate
55   !!----------------------------------------------------------------------   
56
57   SUBROUTINE dyn_hpg( kt )
58      !!---------------------------------------------------------------------
59      !!                  ***  ROUTINE dyn_hpg  ***
60      !!
61      !! ** Purpose :   Compute the now momentum trend due to the hor. gradient
62      !!      of the hydrostatic pressure. Add it to the general momentum trend.
63      !!
64      !! ** Method  :   The now hydrostatic pressure gradient at a given level
65      !!      jk is computed by taking the vertical integral of the in-situ
66      !!      density gradient along the model level from the suface to that
67      !!      level. s-coordinates ('key_s_coord'): a corrective term is added
68      !!      to the horizontal pressure gradient :
69      !!         zhpi = g .....  + 1/e1u mi(rhd) di[ g dep3w ]
70      !!         zhpj = g .....  + 1/e2v mj(rhd) dj[ g dep3w ]
71      !!      add it to the general momentum trend (ua,va).
72      !!         ua = ua - 1/e1u * zhpi
73      !!         va = va - 1/e2v * zhpj
74      !!
75      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
76      !!             - Save the trend in (utrd,vtrd) ('key_trddyn')
77      !!
78      !! History :
79      !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
80      !!        !  91-11  (G. Madec)
81      !!        !  96-01  (G. Madec)  s-coordinates
82      !!        !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg
83      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module, vector opt.
84      !!----------------------------------------------------------------------
85      !! * modules used
86      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
87         &              zhpj => sa      ! use sa as 3D workspace
88
89      !! * Arguments
90      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
91     
92      !! * Local declarations
93      INTEGER ::   ji, jj, jk           ! dummy loop indices
94      REAL(wp) ::   &
95         zcoef0, zcoef1, zuap, zvap     ! temporary scalars
96      !!----------------------------------------------------------------------
97
98      IF( kt == nit000 ) THEN
99         IF(lwp) WRITE(numout,*)
100         IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend'
101         IF(lwp) WRITE(numout,*) '~~~~~~~   s-coordinate case, vector opt. case'
102      ENDIF
103
104      ! 0. Local constant initialization
105      ! --------------------------------
106      zcoef0 = -g * 0.5
107      zuap   = 0.e0
108      zvap   = 0.e0
109
110      ! 1. Surface value
111      ! ----------------
112      DO jj = 2, jpjm1
113         DO ji = fs_2, fs_jpim1   ! vector opt.
114            ! hydrostatic pressure gradient along s-surfaces
115            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj)   &
116                       * ( fse3w(ji+1,jj,1) * rhd(ji+1,jj,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  )
117            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj)   &
118                       * ( fse3w(ji,jj+1,1) * rhd(ji,jj+1,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  )
119            ! s-coordinate pressure gradient correction
120            zuap = -zcoef0 * ( rhd(ji+1,jj,1) + rhd(ji,jj,1) )   &
121                 * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
122            zvap = -zcoef0 * ( rhd(ji,jj+1,1) + rhd(ji,jj,1) )   &
123                 * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
124            ! add to the general momentum trend
125            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
126            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
127# if defined key_trddyn
128            ! save the trend for diagnostics
129            utrd(ji,jj,1,1) = zhpi(ji,jj,1) + zuap
130            vtrd(ji,jj,1,1) = zhpj(ji,jj,1) + zvap
131# endif
132         END DO 
133      END DO 
134
135      ! 2. interior value (2=<jk=<jpkm1)
136      ! -----------------
137      DO jk = 2, jpkm1
138         DO jj = 2, jpjm1 
139            DO ji = fs_2, fs_jpim1   ! vector opt.
140               ! hydrostatic pressure gradient along s-surfaces
141               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
142                           * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
143                              -fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  )
144               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
145                           * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
146                              -fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  )
147               ! s-coordinate pressure gradient correction
148               zuap = -zcoef0 * ( rhd(ji+1,jj  ,jk) + rhd(ji,jj,jk) )   &
149                    * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
150               zvap = -zcoef0 * ( rhd(ji  ,jj+1,jk) + rhd(ji,jj,jk) )   &
151                    * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
152               ! add to the general momentum trend
153               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
154               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
155# if defined key_trddyn
156               ! save the momentum trends for diagnostics
157               utrd(ji,jj,jk,1) = zhpi(ji,jj,jk) + zuap
158               vtrd(ji,jj,jk,1) = zhpj(ji,jj,jk) + zvap
159# endif
160            END DO
161         END DO
162      END DO
163
164      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
165         zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
166         zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
167         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
168         u_ctl = zuap   ;   v_ctl = zvap
169      ENDIF
170
171   END SUBROUTINE dyn_hpg
172
173#elif defined key_partial_steps
174   !!---------------------------------------------------------------------
175   !!   'key_partial_steps'                     z-coordinate partial steps
176   !!---------------------------------------------------------------------
177
178   SUBROUTINE dyn_hpg( kt )
179      !!---------------------------------------------------------------------
180      !!                 ***  ROUTINE dyn_hpg  ***
181      !!                   
182      !! ** Purpose :   Compute the now momentum trend due to the horizontal
183      !!      gradient of the hydrostatic pressure. Add it to the general
184      !!      momentum trend.
185      !!
186      !! ** Method  :   The now hydrostatic pressure gradient at a given level
187      !!      jk is computed by taking the vertical integral of the in-situ
188      !!      density gradient along the model level from the suface to that
189      !!      level:   zhpi = g .....
190      !!               zhpj = g .....
191      !!      add it to the general momentum trend (ua,va).
192      !!            ua = ua - 1/e1u * zhpi
193      !!            va = va - 1/e2v * zhpj
194      !!
195      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
196      !!              - Save the trend in (utrd,vtrd) ('key_trddyn')
197      !!
198      !! History :
199      !!   8.5  !  02-08  (A. Bozec)  Original code
200      !!----------------------------------------------------------------------
201      !! * modules used
202      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
203         &              zhpj => sa      ! use sa as 3D workspace
204
205      !! * Arguments
206      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
207
208      !! * local declarations
209      INTEGER ::   ji, jj, jk           ! dummy loop indices
210      INTEGER ::   iku, ikv             ! temporary integers
211      REAL(wp) ::   &
212         zcoef0, zcoef1, zuap,       &  ! temporary scalars
213         zcoef2, zcoef3, zvap           !    "         "
214      !!----------------------------------------------------------------------
215
216      IF( kt == nit000 ) THEN
217         IF(lwp) WRITE(numout,*)
218         IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend'
219         IF(lwp) WRITE(numout,*) '~~~~~~~   z-coordinate with partial steps'
220         IF(lwp) WRITE(numout,*) '          vector optimization, no autotasking'
221      ENDIF
222
223      ! 0. Local constant initialization
224      ! --------------------------------
225      zcoef0 = -g * 0.5
226      zuap   = 0.e0
227      zvap   = 0.e0
228
229      ! 1. Surface value
230      ! ----------------
231      zcoef1 = zcoef0 * fse3w(1,1,1)
232      DO jj = 2, jpjm1
233         DO ji = fs_2, fs_jpim1   ! vector opt.
234            ! hydrostatic pressure gradient
235            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
236            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
237            ! add to the general momentum trend
238            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
239            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
240# if defined key_trddyn
241            ! save the momentum trends for diagnostics
242            utrd(ji,jj,1,1) = zhpi(ji,jj,1)
243            vtrd(ji,jj,1,1) = zhpj(ji,jj,1)
244# endif
245         END DO
246      END DO
247
248      ! 2. interior value (2=<jk=<jpkm1)
249      ! -----------------
250      DO jk = 2, jpkm1
251         zcoef1 = zcoef0 * fse3w(1,1,jk)
252         DO jj = 2, jpjm1
253            DO ji = fs_2, fs_jpim1   ! vector opt.
254               ! hydrostatic pressure gradient
255               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
256                           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
257                                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
258
259               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
260                           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
261                                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
262               ! add to the general momentum trend
263               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
264               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
265# if defined key_trddyn
266               ! save the momentum trends for diagnostics
267               utrd(ji,jj,jk,1) = zhpi(ji,jj,jk)
268               vtrd(ji,jj,jk,1) = zhpj(ji,jj,jk)
269# endif
270            END DO
271         END DO
272      END DO
273
274      ! partial steps correction at the last level  (new gradient with  intgrd.F)
275# if defined key_vectopt_loop
276         jj = 1
277         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling)
278# else
279      DO jj = 2, jpjm1
280         DO ji = 2, jpim1
281# endif
282            iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1
283            ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
284            zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj  ,iku) )
285            zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji  ,jj+1,ikv) )
286            ! on i-direction
287            IF ( iku > 2 ) THEN
288               ! subtract old value 
289               ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)
290               ! compute the new one   
291               zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1)   &
292                  + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
293               ! add the new one to the general momentum trend
294               ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)
295# if defined key_trddyn
296               ! save the momentum trends for diagnostics
297               utrd(ji,jj,iku,1) = zhpi(ji,jj,iku)
298# endif
299            ENDIF
300            ! on j-direction
301            IF ( ikv > 2 ) THEN
302               ! subtract old value 
303               va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)
304               ! compute the new one   
305               zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1)   &
306                  + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
307               ! add the new one to the general momentum trend
308               va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)
309# if defined key_trddyn
310               ! save the momentum trends for diagnostics
311               vtrd(ji,jj,ikv,1) = zhpj(ji,jj,ikv)
312# endif
313            ENDIF
314# if ! defined key_vectopt_loop
315         END DO
316# endif
317      END DO
318
319      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
320         zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
321         zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
322         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
323         u_ctl = zuap   ;   v_ctl = zvap
324      ENDIF
325
326   END SUBROUTINE dyn_hpg
327
328#else
329   !!---------------------------------------------------------------------
330   !!   Default case :                                        z-coordinate
331   !!---------------------------------------------------------------------
332
333   SUBROUTINE dyn_hpg( kt )
334      !!---------------------------------------------------------------------
335      !!                  ***  ROUTINE dyn_hpg  ***
336      !!
337      !! ** Purpose :   Compute the now momentum trend due to the horizontal
338      !!      gradient of the hydrostatic pressure. Add it to the general
339      !!      momentum trend.
340      !!
341      !! ** Method  :   The now hydrostatic pressure gradient at a given level
342      !!      jk is computed by taking the vertical integral of the in-situ
343      !!      density gradient along the model level from the suface to that
344      !!      level:    zhpi = g .....
345      !!                zhpj = g .....
346      !!      add it to the general momentum trend (ua,va).
347      !!            ua = ua - 1/e1u * zhpi
348      !!            va = va - 1/e2v * zhpj
349      !!
350      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
351      !!             - Save the trend in (utrd,vtrd) ('key_trddyn')
352      !!
353      !! History :
354      !!   1.0  !  87-09  (P. Andrich, m.-a. Foujols)  Original code
355      !!        !  91-11  (G. Madec)
356      !!        !  96-01  (G. Madec)  s-coordinates
357      !!        !  97-05  (G. Madec)  split dynber into dynkeg and dynhpg
358      !!   8.5  !  02-07  (G. Madec)  F90: Free form and module
359      !!----------------------------------------------------------------------
360      !! * modules used
361      USE oce, ONLY :   zhpi => ta,  &  ! use ta as 3D workspace
362         &              zhpj => sa      ! use sa as 3D workspace
363
364      !! * Arguments
365      INTEGER, INTENT( in ) ::   kt     ! ocean time-step index
366
367      !! * local declarations
368      INTEGER ::   ji, jj, jk           ! dummy loop indices
369      REAL(wp) ::   &
370         zcoef0, zcoef1, zuap, zvap     ! temporary scalars
371      !!----------------------------------------------------------------------
372
373      IF( kt == nit000 ) THEN
374         IF(lwp) WRITE(numout,*)
375         IF(lwp) WRITE(numout,*) 'dyn_hpg : hydrostatic pressure gradient trend'
376         IF(lwp) WRITE(numout,*) '~~~~~~~   z-coordinate case '
377      ENDIF
378
379      ! 0. Local constant initialization
380      ! --------------------------------
381      zcoef0 = -g * 0.5
382      zuap   = 0.e0
383      zvap   = 0.e0
384
385      ! 1. Surface value
386      ! ----------------
387      zcoef1 = zcoef0 * fse3w(1,1,1)
388      DO jj = 2, jpjm1
389         DO ji = fs_2, fs_jpim1   ! vector opt.
390            ! hydrostatic pressure gradient
391            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
392            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
393            ! add to the general momentum trend
394            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
395            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
396# if defined key_trddyn
397            ! save the momentum trends for diagnostics
398            utrd(ji,jj,1,1) = zhpi(ji,jj,1)
399            vtrd(ji,jj,1,1) = zhpj(ji,jj,1)
400# endif
401         END DO
402      END DO
403
404      ! 2. interior value (2=<jk=<jpkm1)
405      ! -----------------
406      DO jk = 2, jpkm1
407         zcoef1 = zcoef0 * fse3w(1,1,jk)
408         DO jj = 2, jpjm1
409            DO ji = fs_2, fs_jpim1   ! vector opt.
410               ! hydrostatic pressure gradient
411               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
412                           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
413                                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
414
415               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
416                           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
417                                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
418               ! add to the general momentum trend
419               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
420               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
421# if defined key_trddyn
422               ! save the momentum trends for diagnostics
423               utrd(ji,jj,jk,1) = zhpi(ji,jj,jk)
424               vtrd(ji,jj,jk,1) = zhpj(ji,jj,jk)
425# endif
426            END DO
427         END DO
428      END DO
429
430      IF( l_ctl .AND. lwp ) THEN         ! print sum trends (used for debugging)
431         zuap = SUM( ua(2:jpim1,2:jpjm1,1:jpkm1) * umask(2:jpim1,2:jpjm1,1:jpkm1) )
432         zvap = SUM( va(2:jpim1,2:jpjm1,1:jpkm1) * vmask(2:jpim1,2:jpjm1,1:jpkm1) )
433         WRITE(numout,*) ' hpg  - Ua: ', zuap-u_ctl, ' Va: ', zvap-v_ctl
434         u_ctl = zuap   ;   v_ctl = zvap
435      ENDIF
436
437   END SUBROUTINE dyn_hpg
438
439#endif
440
441   !!======================================================================
442END MODULE dynhpg
Note: See TracBrowser for help on using the repository browser.