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 @ 314

Last change on this file since 314 was 258, checked in by opalod, 19 years ago

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