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

source: trunk/NEMO/OPA_SRC/DYN/dynhpg_jki.F90 @ 466

Last change on this file since 466 was 456, checked in by opalod, 18 years ago

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