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

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.3 KB
Line 
1MODULE dynhpg_jki
2   !!======================================================================
3   !!                       ***  MODULE  dynhpg_jki  ***
4   !! Ocean dynamics:  hydrostatic pressure gradient trend
5   !!======================================================================
6   !! History :  9.0  !  06-09  (G. Madec) From dynhpg module
7   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   hpg_sco_jki : update the momentum trend with the horizontal
11   !!                 gradient of the hydrostatic pressure (s-coordinate)
12   !!   hpg_zps_jki : update the momentum trend with the horizontal
13   !!                 gradient of the hydrostatic pressure (partial step)
14   !!   hpg_zco_jki : update the momentum trend with the horizontal
15   !!                 gradient of the hydrostatic pressure (z-coordinate)
16   !!----------------------------------------------------------------------
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 lbclnk          ! lateral boundary condition
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC   hpg_sco_jki   ! routine called by step.F90
27   PUBLIC   hpg_zps_jki   ! routine called by step.F90
28   PUBLIC   hpg_zco_jki   ! 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   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
37   !!----------------------------------------------------------------------
38
39CONTAINS
40
41   SUBROUTINE hpg_sco_jki( kt )
42      !!---------------------------------------------------------------------
43      !!                  ***  ROUTINE hpg_sco_jki  ***
44      !!       
45      !! ** Purpose :   Compute the now momentum trend due to the horizontal
46      !!     gradient of the hydrostatic pressure. Add it to the general
47      !!     momentum trend.
48      !!
49      !! ** Method  :   The now hydrostatic pressure gradient at a given level
50      !!      jk is computed by taking the vertical integral of the in-situ
51      !!      density gradient along the model level from the suface to that
52      !!      level. s-coordinate case (ln_sco=T): a corrective term is
53      !!      added to the horizontal pressure gradient :
54      !!         zhpi = grav .....   + 1/e1u mi(rhd) di[ grav dep3w ]
55      !!         zhpj = grav .....   + 1/e2v mj(rhd) dj[ grav dep3w ]
56      !!      add it to the general momentum trend (ua,va).
57      !!         ua = ua - 1/e1u * zhpi
58      !!         va = va - 1/e2v * zhpj
59      !!      j-k-i loop (j-slab) ('key_mpp_omp')
60      !!
61      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
62      !!             - Save the trend in (ztrdu,ztrdv) ('key_trddyn')
63      !!----------------------------------------------------------------------
64      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
65      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
66      !!
67      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
68      !!
69      INTEGER  ::   ji, jj, jk           ! dummy loop indices
70      REAL(wp) ::   zcoef0, zuap, zvap   ! temporary scalars
71      !!----------------------------------------------------------------------
72
73      IF( kt == nit000 ) THEN
74         IF(lwp) WRITE(numout,*)
75         IF(lwp) WRITE(numout,*) 'hpg_sco_jki : s-coordinate hydrostatic pressure gradient trend'
76         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   OpenMP / NEC autotasking case (j-k-i loop)'
77      ENDIF
78
79      ! Local constant initialization
80      zcoef0 = - grav * 0.5
81      zuap   = 0.e0
82      zvap   = 0.e0
83
84      !                                                ! ===============
85      DO jj = 2, jpjm1                                 !  Vertical slab
86         !                                             ! ===============
87         DO ji = 2, jpim1         !  Surface value
88            ! hydrostatic pressure gradient along s-surfaces
89            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj)   &
90                       * ( fse3w(ji+1,jj,1) * rhd(ji+1,jj,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  )
91            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj)   &
92                       * ( fse3w(ji,jj+1,1) * rhd(ji,jj+1,1) - fse3w(ji,jj,1) * rhd(ji,jj,1)  )
93            ! s-coordinate pressure gradient correction
94            zuap = -zcoef0 * ( rhd(ji+1,jj,1) + rhd(ji,jj,1) )   &
95                 * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
96            zvap = -zcoef0 * ( rhd(ji,jj+1,1) + rhd(ji,jj,1) )   &
97                 * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
98            ! add to the general momentum trend
99            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
100            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
101         END DO 
102         !
103         DO jk = 2, jpkm1         ! interior value (2=<jk=<jpkm1)
104            DO ji = 2, jpim1
105               ! hydrostatic pressure gradient along s-surfaces
106               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
107                  &           * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
108                  &              -fse3w(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  )
109               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
110                  &           * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
111                  &              -fse3w(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  )
112               ! s-coordinate pressure gradient correction
113               zuap = -zcoef0 * ( rhd(ji+1,jj  ,jk) + rhd(ji,jj,jk) )   &
114                    * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
115               zvap = -zcoef0 * ( rhd(ji  ,jj+1,jk) + rhd(ji,jj,jk) )   &
116                    * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
117               ! add to the general momentum trend
118               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
119               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
120            END DO
121         END DO
122         !                                             ! ===============
123      END DO                                           !   End of slab
124      !                                                ! ===============
125   END SUBROUTINE hpg_sco_jki
126
127
128   SUBROUTINE hpg_zps_jki( kt )
129      !!---------------------------------------------------------------------
130      !!                  ***  ROUTINE hpg_zps_jki  ***
131      !! 
132      !! ** Purpose :   Compute the now momentum trend due to the hor. gradient
133      !!      of the hydrostatic pressure. Add it to the general momentum trend.
134      !!
135      !! ** Method  :   The now hydrostatic pressure gradient at a given level
136      !!      jk is computed by taking the vertical integral of the in-situ
137      !!      density gradient along the model level from the suface to that
138      !!      level:    zhpi = grav .....
139      !!                zhpj = grav .....
140      !!      add it to the general momentum trend (ua,va).
141      !!            ua = ua - 1/e1u * zhpi
142      !!            va = va - 1/e2v * zhpj
143      !!      j-k-i loop (j-slab) ('key_mpp_omp')
144      !!
145      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
146      !!             - Save the trend in (ztrdu,ztrdv) ('key_trddyn')
147      !!----------------------------------------------------------------------
148      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
149      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
150      !!
151      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
152      !!
153      INTEGER  ::   ji, jj, jk             ! dummy loop indices
154      INTEGER  ::   iku, ikv               ! temporary integers
155      REAL(wp) ::   zcoef0, zcoef1, zuap   ! temporary scalars
156      REAL(wp) ::   zcoef2, zcoef3, zvap   !    "         "
157      !!----------------------------------------------------------------------
158
159      IF( kt == nit000 ) THEN
160         IF(lwp) WRITE(numout,*)
161         IF(lwp) WRITE(numout,*) 'hpg_zps_jki : z-coord. partial steps hydrostatic pressure gradient trend'
162         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   OpenMP / NEC autotasking case (j-k-i loop)'
163      ENDIF
164
165      ! Local constant initialization
166      zcoef0 = - grav * 0.5
167      zuap   = 0.e0
168      zvap   = 0.e0
169      !                                                ! ===============
170      DO jj = 2, jpjm1                                 !  Vertical slab
171         !                                             ! ===============
172         DO ji = 2, jpim1         ! Surface value
173            zcoef1 = zcoef0 * fse3w(ji,jj,1)
174            ! hydrostatic pressure gradient
175            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
176            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
177            ! add to the general momentum trend
178            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
179            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
180         END DO
181         !
182         DO jk = 2, jpkm1         ! interior value (2=<jk=<jpkm1)
183            DO ji = 2, jpim1
184               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
185               ! hydrostatic pressure gradient
186               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
187                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
188                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
189
190               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
191                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
192                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
193               ! add to the general momentum trend
194               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
195               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
196            END DO
197         END DO
198         !
199         ! partial steps correction at the last level  (new gradient with  intgrd.F)
200         DO ji = 2, jpim1
201            iku = MIN ( mbathy(ji,jj), mbathy(ji+1,jj) ) - 1
202            ikv = MIN ( mbathy(ji,jj), mbathy(ji,jj+1) ) - 1
203            zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj  ,iku) )
204            zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji  ,jj+1,ikv) )
205            ! on i-direction
206            IF ( iku > 2 ) THEN
207               ! subtract old value 
208               ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)
209               ! compute the new one   
210               zhpi (ji,jj,iku) = zhpi(ji,jj,iku-1)   &
211                  + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
212               ! add the new one to the general momentum trend
213               ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)
214            ENDIF
215            ! on j-direction
216            IF ( ikv > 2 ) THEN
217               ! subtract old value 
218               va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)
219               ! compute the new one   
220               zhpj (ji,jj,ikv) = zhpj(ji,jj,ikv-1)   &
221                  + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
222               ! add the new one to the general momentum trend
223               va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)
224            ENDIF
225         END DO
226         !                                             ! ===============
227      END DO                                           !   End of slab
228      !                                                ! ===============
229   END SUBROUTINE hpg_zps_jki
230
231
232   SUBROUTINE hpg_zco_jki( kt )
233      !!---------------------------------------------------------------------
234      !!                  ***  ROUTINE hpg_zco_jki  ***
235      !!   
236      !! ** Purpose :   Compute the now momentum trend due to the horizontal
237      !!      gradient of the hydrostatic pressure. Add it to the general
238      !!      momentum trend.
239      !!
240      !! ** Method  :   The now hydrostatic pressure gradient at a given level
241      !!      jk is computed by taking the vertical integral of the in-situ 
242      !!      density gradient along the model level from the suface to that
243      !!      level:    zhpi = grav .....
244      !!                zhpj = grav .....
245      !!      add it to the general momentum trend (ua,va).
246      !!            ua = ua - 1/e1u * zhpi
247      !!            va = va - 1/e2v * zhpj
248      !!      j-k-i loop (j-slab) ('key_mpp_omp')
249      !!
250      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
251      !!             - Save the trend in (ztrdu,ztrdv) ('key_trddyn')
252      !!----------------------------------------------------------------------
253      USE oce, ONLY :   zhpi => ta   ! use ta as 3D workspace
254      USE oce, ONLY :   zhpj => sa   ! use sa as 3D workspace
255      !!
256      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
257      !!
258      INTEGER  ::   ji, jj, jk                   ! dummy loop indices
259      REAL(wp) ::   zcoef0, zcoef1, zuap, zvap   ! temporary scalars
260      !!----------------------------------------------------------------------
261
262      IF( kt == nit000 ) THEN
263         IF(lwp) WRITE(numout,*)
264         IF(lwp) WRITE(numout,*) 'hpg_zco_jki : z-coordinate hydrostatic pressure gradient trend'
265         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   OpenMP / NEC auto-tasking (j-k-i loop)'
266      ENDIF
267
268      ! Local constant initialization
269      zcoef0 = - grav * 0.5
270      zuap   = 0.e0
271      zvap   = 0.e0
272
273      !                                                ! ===============
274      DO jj = 2, jpjm1                                 !  Vertical slab
275         !                                             ! ===============
276         DO ji = 2, jpim1         ! Surface value
277            zcoef1 = zcoef0 * fse3w(ji,jj,1)
278            ! hydrostatic pressure gradient
279            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
280            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
281            ! add to the general momentum trend
282            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
283            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
284         END DO
285         !
286         DO jk = 2, jpkm1         ! interior value (2=<jk=<jpkm1)
287            DO ji = 2, jpim1
288               zcoef1 = zcoef0 * fse3w(ji,jj,jk)
289               ! hydrostatic pressure gradient
290               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
291                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )   &
292                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) / e1u(ji,jj)
293
294               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
295                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )   &
296                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) / e2v(ji,jj)
297               ! add to the general momentum trend
298               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
299               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
300            END DO
301         END DO
302         !                                             ! ===============
303      END DO                                           !   End of slab
304      !                                                ! ===============
305   END SUBROUTINE hpg_zco_jki
306
307   !!======================================================================
308END MODULE dynhpg_jki
Note: See TracBrowser for help on using the repository browser.