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

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

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