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

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

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