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

source: trunk/NEMO/OPA_SRC/DYN/dynhpg_atsk.F90 @ 78

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

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