source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

  • Property svn:keywords set to Id
File size: 70.0 KB
Line 
1MODULE dynhpg
2   !!======================================================================
3   !!                       ***  MODULE  dynhpg  ***
4   !! Ocean dynamics:  hydrostatic pressure gradient trend
5   !!======================================================================
6   !! History :  OPA  !  1987-09  (P. Andrich, M.-A. Foujols)  hpg_zco: Original code
7   !!            5.0  !  1991-11  (G. Madec)
8   !!            7.0  !  1996-01  (G. Madec)  hpg_sco: Original code for s-coordinates
9   !!            8.0  !  1997-05  (G. Madec)  split dynber into dynkeg and dynhpg
10   !!            8.5  !  2002-07  (G. Madec)  F90: Free form and module
11   !!            8.5  !  2002-08  (A. Bozec)  hpg_zps: Original code
12   !!   NEMO     1.0  !  2005-10  (A. Beckmann, B.W. An)  various s-coordinate options
13   !!                 !         Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot
14   !!             -   !  2005-11  (G. Madec) style & small optimisation
15   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
16   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates
17   !!                 !           (A. Coward) suppression of hel, wdj and rot options
18   !!            3.6  !  2014-11  (P. Mathiot) hpg_isf: original code for ice shelf cavity
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dyn_hpg      : update the momentum trend with the now horizontal
23   !!                  gradient of the hydrostatic pressure
24   !!   dyn_hpg_init : initialisation and control of options
25   !!       hpg_zco  : z-coordinate scheme
26   !!       hpg_zps  : z-coordinate plus partial steps (interpolation)
27   !!       hpg_sco  : s-coordinate (standard jacobian formulation)
28   !!       hpg_isf  : s-coordinate (sco formulation) adapted to ice shelf
29   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial)
30   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial)
31   !!----------------------------------------------------------------------
32   USE oce             ! ocean dynamics and tracers
33   USE sbc_oce         ! surface variable (only for the flag with ice shelf)
34   USE dom_oce         ! ocean space and time domain
35   USE wet_dry         ! wetting and drying
36   USE phycst          ! physical constants
37   USE trd_oce         ! trends: ocean variables
38   USE trddyn          ! trend manager: dynamics
39!jc   USE zpshde          ! partial step: hor. derivative     (zps_hde routine)
40   !
41   USE in_out_manager  ! I/O manager
42   USE prtctl          ! Print control
43   USE lbclnk          ! lateral boundary condition
44   USE lib_mpp         ! MPP library
45   USE eosbn2          ! compute density
46   USE timing          ! Timing
47   USE iom
48
49   IMPLICIT NONE
50   PRIVATE
51
52   PUBLIC   dyn_hpg        ! routine called by step module
53   PUBLIC   dyn_hpg_init   ! routine called by opa module
54
55   !                                 !!* Namelist namdyn_hpg : hydrostatic pressure gradient
56   LOGICAL , PUBLIC ::   ln_hpg_zco   !: z-coordinate - full steps
57   LOGICAL , PUBLIC ::   ln_hpg_zps   !: z-coordinate - partial steps (interpolation)
58   LOGICAL , PUBLIC ::   ln_hpg_sco   !: s-coordinate (standard jacobian formulation)
59   LOGICAL , PUBLIC ::   ln_hpg_djc   !: s-coordinate (Density Jacobian with Cubic polynomial)
60   LOGICAL , PUBLIC ::   ln_hpg_prj   !: s-coordinate (Pressure Jacobian scheme)
61   LOGICAL , PUBLIC ::   ln_hpg_isf   !: s-coordinate similar to sco modify for isf
62
63   INTEGER , PUBLIC ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
64
65   !! * Substitutions
66#  include "vectopt_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
69   !! $Id$
70   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
71   !!----------------------------------------------------------------------
72CONTAINS
73
74   SUBROUTINE dyn_hpg( kt )
75      !!---------------------------------------------------------------------
76      !!                  ***  ROUTINE dyn_hpg  ***
77      !!
78      !! ** Method  :   Call the hydrostatic pressure gradient routine
79      !!              using the scheme defined in the namelist
80      !!
81      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
82      !!             - send trends to trd_dyn for futher diagnostics (l_trddyn=T)
83      !!----------------------------------------------------------------------
84      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
85      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztrdu, ztrdv
86      !!----------------------------------------------------------------------
87      !
88      IF( nn_timing == 1 )  CALL timing_start('dyn_hpg')
89      !
90      IF( l_trddyn ) THEN                    ! Temporary saving of ua and va trends (l_trddyn)
91         ztrdu(:,:,:) = ua(:,:,:)
92         ztrdv(:,:,:) = va(:,:,:)
93      ENDIF
94      !
95      SELECT CASE ( nhpg )      ! Hydrostatic pressure gradient computation
96      CASE (  0 )   ;   CALL hpg_zco    ( kt )      ! z-coordinate
97      CASE (  1 )   ;   CALL hpg_zps    ( kt )      ! z-coordinate plus partial steps (interpolation)
98      CASE (  2 )   ;   CALL hpg_sco    ( kt )      ! s-coordinate (standard jacobian formulation)
99      CASE (  3 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial)
100      CASE (  4 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme)
101      CASE (  5 )   ;   CALL hpg_isf    ( kt )      ! s-coordinate similar to sco modify for ice shelf
102      END SELECT
103      !
104      IF( l_trddyn ) THEN      ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
105         ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
106         ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
107         CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )
108      ENDIF
109      !
110      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg  - Ua: ', mask1=umask,   &
111         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
112      !
113      IF( nn_timing == 1 )  CALL timing_stop('dyn_hpg')
114      !
115   END SUBROUTINE dyn_hpg
116
117
118   SUBROUTINE dyn_hpg_init
119      !!----------------------------------------------------------------------
120      !!                 ***  ROUTINE dyn_hpg_init  ***
121      !!
122      !! ** Purpose :   initializations for the hydrostatic pressure gradient
123      !!              computation and consistency control
124      !!
125      !! ** Action  :   Read the namelist namdyn_hpg and check the consistency
126      !!      with the type of vertical coordinate used (zco, zps, sco)
127      !!----------------------------------------------------------------------
128      INTEGER ::   ioptio = 0      ! temporary integer
129      INTEGER ::   ios             ! Local integer output status for namelist read
130      !!
131      INTEGER  ::   ji, jj, jk, ikt    ! dummy loop indices      ISF
132      REAL(wp) ::   znad
133      REAL(wp), DIMENSION(jpi,jpj,2)   ::  ztstop       ! hypothesys on isf density
134      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::  zrhd         ! hypothesys on isf density
135      REAL(wp), DIMENSION(jpi,jpj)     ::  zrhdtop_isf  ! density at bottom of ISF
136      REAL(wp), DIMENSION(jpi,jpj)     ::  ziceload     ! density at bottom of ISF
137      !!
138      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco,     &
139         &                 ln_hpg_djc, ln_hpg_prj, ln_hpg_isf
140      !!----------------------------------------------------------------------
141      !
142      REWIND( numnam_ref )              ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient
143      READ  ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901)
144901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp )
145      !
146      REWIND( numnam_cfg )              ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient
147      READ  ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 )
148902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )
149      IF(lwm) WRITE ( numond, namdyn_hpg )
150      !
151      IF(lwp) THEN                   ! Control print
152         WRITE(numout,*)
153         WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
154         WRITE(numout,*) '~~~~~~~~~~~~'
155         WRITE(numout,*) '   Namelist namdyn_hpg : choice of hpg scheme'
156         WRITE(numout,*) '      z-coord. - full steps                             ln_hpg_zco    = ', ln_hpg_zco
157         WRITE(numout,*) '      z-coord. - partial steps (interpolation)          ln_hpg_zps    = ', ln_hpg_zps
158         WRITE(numout,*) '      s-coord. (standard jacobian formulation)          ln_hpg_sco    = ', ln_hpg_sco
159         WRITE(numout,*) '      s-coord. (standard jacobian formulation) for isf  ln_hpg_isf    = ', ln_hpg_isf
160         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc
161         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj
162      ENDIF
163      !
164      IF( ln_hpg_djc )   &
165         &   CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method &
166                           & currently disabled (bugs under investigation). Please select &
167                           & either  ln_hpg_sco or  ln_hpg_prj instead')
168      !
169      IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) )        &
170         &   CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', &
171         &                 '   the standard jacobian formulation hpg_sco    or '    , &
172         &                 '   the pressure jacobian formulation hpg_prj'            )
173
174      IF(       ln_hpg_isf .AND. .NOT. ln_isfcav )   &
175         &   CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' )
176      IF( .NOT. ln_hpg_isf .AND.       ln_isfcav )   &
177         &   CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
178      !
179      !                               ! Set nhpg from ln_hpg_... flags
180      IF( ln_hpg_zco )   nhpg = 0
181      IF( ln_hpg_zps )   nhpg = 1
182      IF( ln_hpg_sco )   nhpg = 2
183      IF( ln_hpg_djc )   nhpg = 3
184      IF( ln_hpg_prj )   nhpg = 4
185      IF( ln_hpg_isf )   nhpg = 5
186      !
187      !                               ! Consistency check
188      ioptio = 0
189      IF( ln_hpg_zco )   ioptio = ioptio + 1
190      IF( ln_hpg_zps )   ioptio = ioptio + 1
191      IF( ln_hpg_sco )   ioptio = ioptio + 1
192      IF( ln_hpg_djc )   ioptio = ioptio + 1
193      IF( ln_hpg_prj )   ioptio = ioptio + 1
194      IF( ln_hpg_isf )   ioptio = ioptio + 1
195      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
196      !
197      ! initialisation of ice shelf load
198      IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0
199      IF (       ln_isfcav ) THEN
200         !
201         IF(lwp) WRITE(numout,*)
202         IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
203         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 
204
205         ! To use density and not density anomaly
206         znad=1._wp
207         
208         ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
209         ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
210
211         ! compute density of the water displaced by the ice shelf
212         DO jk = 1, jpk
213            CALL eos(ztstop(:,:,:),gdept_n(:,:,jk),zrhd(:,:,jk))
214         END DO
215     
216         ! compute rhd at the ice/oce interface (ice shelf side)
217         CALL eos(ztstop,risfdep,zrhdtop_isf)
218
219         ! Surface value + ice shelf gradient
220         ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
221         ! divided by 2 later
222         ziceload = 0._wp
223         DO jj = 1, jpj
224            DO ji = 1, jpi
225               ikt=mikt(ji,jj)
226               ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1))
227               DO jk=2,ikt-1
228                  ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * e3w_n(ji,jj,jk) &
229                     &                              * (1._wp - tmask(ji,jj,jk))
230               END DO
231               IF (ikt  >=  2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) &
232                                  &                              * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
233            END DO
234         END DO
235         riceload(:,:)=ziceload(:,:)  ! need to be saved for diaar5
236
237      END IF
238      !
239   END SUBROUTINE dyn_hpg_init
240
241
242   SUBROUTINE hpg_zco( kt )
243      !!---------------------------------------------------------------------
244      !!                  ***  ROUTINE hpg_zco  ***
245      !!
246      !! ** Method  :   z-coordinate case, levels are horizontal surfaces.
247      !!      The now hydrostatic pressure gradient at a given level, jk,
248      !!      is computed by taking the vertical integral of the in-situ
249      !!      density gradient along the model level from the suface to that
250      !!      level:    zhpi = grav .....
251      !!                zhpj = grav .....
252      !!      add it to the general momentum trend (ua,va).
253      !!            ua = ua - 1/e1u * zhpi
254      !!            va = va - 1/e2v * zhpj
255      !!
256      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
257      !!----------------------------------------------------------------------
258      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
259      !
260      INTEGER  ::   ji, jj, jk       ! dummy loop indices
261      REAL(wp) ::   zcoef0, zcoef1   ! temporary scalars
262      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
263      !!----------------------------------------------------------------------
264      !
265      !
266      IF( kt == nit000 ) THEN
267         IF(lwp) WRITE(numout,*)
268         IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
269         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate case '
270      ENDIF
271
272      zcoef0 = - grav * 0.5_wp      ! Local constant initialization
273
274      ! Surface value
275      DO jj = 2, jpjm1
276         DO ji = fs_2, fs_jpim1   ! vector opt.
277            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
278            ! hydrostatic pressure gradient
279            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
280            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_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      END DO
286
287      !
288      ! interior value (2=<jk=<jpkm1)
289      DO jk = 2, jpkm1
290         DO jj = 2, jpjm1
291            DO ji = fs_2, fs_jpim1   ! vector opt.
292               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
293               ! hydrostatic pressure gradient
294               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
295                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) )    &
296                  &                       - ( rhd(ji  ,jj,jk)+rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
297
298               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
299                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) )    &
300                  &                       - ( rhd(ji,jj,  jk)+rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
301               ! add to the general momentum trend
302               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
303               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
304            END DO
305         END DO
306      END DO
307      !
308      !
309   END SUBROUTINE hpg_zco
310
311
312   SUBROUTINE hpg_zps( kt )
313      !!---------------------------------------------------------------------
314      !!                 ***  ROUTINE hpg_zps  ***
315      !!
316      !! ** Method  :   z-coordinate plus partial steps case.  blahblah...
317      !!
318      !! ** Action  : - Update (ua,va) with the now hydrastatic pressure trend
319      !!----------------------------------------------------------------------
320      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
321      !!
322      INTEGER  ::   ji, jj, jk                       ! dummy loop indices
323      INTEGER  ::   iku, ikv                         ! temporary integers
324      REAL(wp) ::   zcoef0, zcoef1, zcoef2, zcoef3   ! temporary scalars
325      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
326      !!----------------------------------------------------------------------
327      !
328      !
329      IF( kt == nit000 ) THEN
330         IF(lwp) WRITE(numout,*)
331         IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
332         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   z-coordinate with partial steps - vector optimization'
333      ENDIF
334
335      ! Partial steps: bottom before horizontal gradient of t, s, rd at the last ocean level
336!jc      CALL zps_hde    ( kt, jpts, tsn, gtsu, gtsv, rhd, gru , grv )
337
338      ! Local constant initialization
339      zcoef0 = - grav * 0.5_wp
340
341      !  Surface value (also valid in partial step case)
342      DO jj = 2, jpjm1
343         DO ji = fs_2, fs_jpim1   ! vector opt.
344            zcoef1 = zcoef0 * e3w_n(ji,jj,1)
345            ! hydrostatic pressure gradient
346            zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj  ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
347            zhpj(ji,jj,1) = zcoef1 * ( rhd(ji  ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
348            ! add to the general momentum trend
349            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
350            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
351         END DO
352      END DO
353
354      ! interior value (2=<jk=<jpkm1)
355      DO jk = 2, jpkm1
356         DO jj = 2, jpjm1
357            DO ji = fs_2, fs_jpim1   ! vector opt.
358               zcoef1 = zcoef0 * e3w_n(ji,jj,jk)
359               ! hydrostatic pressure gradient
360               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)   &
361                  &           + zcoef1 * (  ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) )   &
362                  &                       - ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) )  ) * r1_e1u(ji,jj)
363
364               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)   &
365                  &           + zcoef1 * (  ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) )   &
366                  &                       - ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) )  ) * r1_e2v(ji,jj)
367               ! add to the general momentum trend
368               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
369               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
370            END DO
371         END DO
372      END DO
373
374      ! partial steps correction at the last level  (use gru & grv computed in zpshde.F90)
375      DO jj = 2, jpjm1
376         DO ji = 2, jpim1
377            iku = mbku(ji,jj)
378            ikv = mbkv(ji,jj)
379            zcoef2 = zcoef0 * MIN( e3w_n(ji,jj,iku), e3w_n(ji+1,jj  ,iku) )
380            zcoef3 = zcoef0 * MIN( e3w_n(ji,jj,ikv), e3w_n(ji  ,jj+1,ikv) )
381            IF( iku > 1 ) THEN            ! on i-direction (level 2 or more)
382               ua  (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku)         ! subtract old value
383               zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1)                   &   ! compute the new one
384                  &            + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj)
385               ua  (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku)         ! add the new one to the general momentum trend
386            ENDIF
387            IF( ikv > 1 ) THEN            ! on j-direction (level 2 or more)
388               va  (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv)         ! subtract old value
389               zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1)                   &   ! compute the new one
390                  &            + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj)
391               va  (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv)         ! add the new one to the general momentum trend
392            ENDIF
393         END DO
394      END DO
395      !
396      !
397   END SUBROUTINE hpg_zps
398
399
400   SUBROUTINE hpg_sco( kt )
401      !!---------------------------------------------------------------------
402      !!                  ***  ROUTINE hpg_sco  ***
403      !!
404      !! ** Method  :   s-coordinate case. Jacobian scheme.
405      !!      The now hydrostatic pressure gradient at a given level, jk,
406      !!      is computed by taking the vertical integral of the in-situ
407      !!      density gradient along the model level from the suface to that
408      !!      level. s-coordinates (ln_sco): a corrective term is added
409      !!      to the horizontal pressure gradient :
410      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
411      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
412      !!      add it to the general momentum trend (ua,va).
413      !!         ua = ua - 1/e1u * zhpi
414      !!         va = va - 1/e2v * zhpj
415      !!
416      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
417      !!----------------------------------------------------------------------
418      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
419      !!
420      INTEGER  ::   ji, jj, jk, jii, jjj                 ! dummy loop indices
421      REAL(wp) ::   zcoef0, zuap, zvap, znad, ztmp       ! temporary scalars
422      LOGICAL  ::   ll_tmp1, ll_tmp2                     ! local logical variables
423      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
424      REAL(wp), DIMENSION(jpi,jpj)   ::  zcpx, zcpy !W/D pressure filter
425      !!----------------------------------------------------------------------
426      !
427      !
428      IF( kt == nit000 ) THEN
429         IF(lwp) WRITE(numout,*)
430         IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
431         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, OPA original scheme used'
432      ENDIF
433      !
434      zcoef0 = - grav * 0.5_wp
435      IF ( ln_linssh ) THEN   ;   znad = 0._wp         ! Fixed    volume: density anomaly
436      ELSE                    ;   znad = 1._wp         ! Variable volume: density
437      ENDIF
438      !
439      IF( ln_wd ) THEN
440        DO jj = 2, jpjm1
441           DO ji = 2, jpim1 
442             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
443                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
444                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
445                  &                                                         > rn_wdmin1 + rn_wdmin2
446             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
447                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
448                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
449
450             IF(ll_tmp1) THEN
451               zcpx(ji,jj) = 1.0_wp
452             ELSE IF(ll_tmp2) THEN
453               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
454               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
455                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
456             ELSE
457               zcpx(ji,jj) = 0._wp
458             END IF
459     
460             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
461                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
462                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
463                  &                                                         > rn_wdmin1 + rn_wdmin2
464             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
465                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
466                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
467
468             IF(ll_tmp1) THEN
469               zcpy(ji,jj) = 1.0_wp
470             ELSE IF(ll_tmp2) THEN
471               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
472               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
473                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
474             ELSE
475               zcpy(ji,jj) = 0._wp
476             END IF
477           END DO
478        END DO
479        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
480      END IF
481
482      ! Surface value
483      DO jj = 2, jpjm1
484         DO ji = fs_2, fs_jpim1   ! vector opt.
485            ! hydrostatic pressure gradient along s-surfaces
486            zhpi(ji,jj,1) = zcoef0 * (  e3w_n(ji+1,jj  ,1) * ( znad + rhd(ji+1,jj  ,1) )    &
487               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e1u(ji,jj)
488            zhpj(ji,jj,1) = zcoef0 * (  e3w_n(ji  ,jj+1,1) * ( znad + rhd(ji  ,jj+1,1) )    &
489               &                      - e3w_n(ji  ,jj  ,1) * ( znad + rhd(ji  ,jj  ,1) )  ) * r1_e2v(ji,jj)
490            ! s-coordinate pressure gradient correction
491            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
492               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
493            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
494               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
495
496
497            IF( ln_wd ) THEN
498
499              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
500              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
501              zuap = zuap * zcpx(ji,jj)
502              zvap = zvap * zcpy(ji,jj)
503            ENDIF
504
505            ! add to the general momentum trend
506            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
507            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
508         END DO
509      END DO
510
511      ! interior value (2=<jk=<jpkm1)
512      DO jk = 2, jpkm1
513         DO jj = 2, jpjm1
514            DO ji = fs_2, fs_jpim1   ! vector opt.
515               ! hydrostatic pressure gradient along s-surfaces
516               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj)   &
517                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )   &
518                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad )  )
519               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj)   &
520                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad )   &
521                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad )  )
522               ! s-coordinate pressure gradient correction
523               zuap = -zcoef0 * ( rhd    (ji+1,jj  ,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
524                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj)
525               zvap = -zcoef0 * ( rhd    (ji  ,jj+1,jk) + rhd    (ji,jj,jk) + 2._wp * znad )   &
526                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj)
527
528               IF( ln_wd ) THEN
529                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
530                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
531                 zuap = zuap * zcpx(ji,jj)
532                 zvap = zvap * zcpy(ji,jj)
533               ENDIF
534
535               ! add to the general momentum trend
536               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
537               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
538            END DO
539         END DO
540      END DO
541      !
542      !
543   END SUBROUTINE hpg_sco
544
545
546   SUBROUTINE hpg_isf( kt )
547      !!---------------------------------------------------------------------
548      !!                  ***  ROUTINE hpg_isf  ***
549      !!
550      !! ** Method  :   s-coordinate case. Jacobian scheme.
551      !!      The now hydrostatic pressure gradient at a given level, jk,
552      !!      is computed by taking the vertical integral of the in-situ
553      !!      density gradient along the model level from the suface to that
554      !!      level. s-coordinates (ln_sco): a corrective term is added
555      !!      to the horizontal pressure gradient :
556      !!         zhpi = grav .....  + 1/e1u mi(rhd) di[ grav dep3w ]
557      !!         zhpj = grav .....  + 1/e2v mj(rhd) dj[ grav dep3w ]
558      !!      add it to the general momentum trend (ua,va).
559      !!         ua = ua - 1/e1u * zhpi
560      !!         va = va - 1/e2v * zhpj
561      !!      iceload is added and partial cell case are added to the top and bottom
562      !!     
563      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
564      !!----------------------------------------------------------------------
565      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
566      !!
567      INTEGER  ::   ji, jj, jk, ikt, iktp1i, iktp1j   ! dummy loop indices
568      REAL(wp) ::   zcoef0, zuap, zvap, znad          ! temporary scalars
569      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::  zhpi, zhpj
570      REAL(wp), DIMENSION(jpi,jpj,2)   ::  ztstop
571      REAL(wp), DIMENSION(jpi,jpj)     ::  zrhdtop_oce
572      !!----------------------------------------------------------------------
573      !
574      !
575      ! Local constant initialization
576      zcoef0 = - grav * 0.5_wp
577 
578      ! To use density and not density anomaly
579      znad=1._wp
580
581      ! iniitialised to 0. zhpi zhpi
582      zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
583
584      ! compute rhd at the ice/oce interface (ocean side)
585      ! usefull to reduce residual current in the test case ISOMIP with no melting
586      DO ji=1,jpi
587        DO jj=1,jpj
588          ikt=mikt(ji,jj)
589          ztstop(ji,jj,1)=tsn(ji,jj,ikt,1)
590          ztstop(ji,jj,2)=tsn(ji,jj,ikt,2)
591        END DO
592      END DO
593      CALL eos( ztstop, risfdep, zrhdtop_oce )
594
595!==================================================================================     
596!===== Compute surface value =====================================================
597!==================================================================================
598      DO jj = 2, jpjm1
599         DO ji = fs_2, fs_jpim1   ! vector opt.
600            ikt    = mikt(ji,jj)
601            iktp1i = mikt(ji+1,jj)
602            iktp1j = mikt(ji,jj+1)
603            ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
604            ! we assume ISF is in isostatic equilibrium
605            zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i)                                    &
606               &                                    * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) )   &
607               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         &
608               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
609               &                                  + ( riceload(ji+1,jj) - riceload(ji,jj))                            ) 
610            zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j)                                    &
611               &                                    * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) )   &
612               &                                  - 0.5_wp * e3w_n(ji,jj,ikt)                                         & 
613               &                                    * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) )          &
614               &                                  + ( riceload(ji,jj+1) - riceload(ji,jj))                            ) 
615            ! s-coordinate pressure gradient correction (=0 if z coordinate)
616            zuap = -zcoef0 * ( rhd    (ji+1,jj,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
617               &           * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj)
618            zvap = -zcoef0 * ( rhd    (ji,jj+1,1) + rhd    (ji,jj,1) + 2._wp * znad )   &
619               &           * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj)
620            ! add to the general momentum trend
621            ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
622            va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
623         END DO
624      END DO
625!==================================================================================     
626!===== Compute interior value =====================================================
627!==================================================================================
628      ! interior value (2=<jk=<jpkm1)
629      DO jk = 2, jpkm1
630         DO jj = 2, jpjm1
631            DO ji = fs_2, fs_jpim1   ! vector opt.
632               ! hydrostatic pressure gradient along s-surfaces
633               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj)   &
634                  &           * (  e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk)   &
635                  &              - e3w_n(ji  ,jj,jk) * ( rhd(ji  ,jj,jk) + rhd(ji  ,jj,jk-1) + 2*znad ) * wmask(ji  ,jj,jk)   )
636               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj)   &
637                  &           * (  e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk)   &
638                  &              - e3w_n(ji,jj  ,jk) * ( rhd(ji,jj,  jk) + rhd(ji,jj  ,jk-1) + 2*znad ) * wmask(ji,jj  ,jk)   )
639               ! s-coordinate pressure gradient correction
640               zuap = -zcoef0 * ( rhd   (ji+1,jj  ,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
641                  &           * ( gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj)
642               zvap = -zcoef0 * ( rhd   (ji  ,jj+1,jk) + rhd   (ji,jj,jk) + 2._wp * znad )   &
643                  &           * ( gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj)
644               ! add to the general momentum trend
645               ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
646               va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
647            END DO
648         END DO
649      END DO
650     !
651      !
652   END SUBROUTINE hpg_isf
653
654
655   SUBROUTINE hpg_djc( kt )
656      !!---------------------------------------------------------------------
657      !!                  ***  ROUTINE hpg_djc  ***
658      !!
659      !! ** Method  :   Density Jacobian with Cubic polynomial scheme
660      !!
661      !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
662      !!----------------------------------------------------------------------
663      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
664      !!
665      INTEGER  ::   ji, jj, jk          ! dummy loop indices
666      REAL(wp) ::   zcoef0, zep, cffw   ! temporary scalars
667      REAL(wp) ::   z1_10, cffu, cffx   !    "         "
668      REAL(wp) ::   z1_12, cffv, cffy   !    "         "
669      LOGICAL  ::   ll_tmp1, ll_tmp2    ! local logical variables
670      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zhpi, zhpj
671      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  dzx, dzy, dzz, dzu, dzv, dzw
672      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  drhox, drhoy, drhoz, drhou, drhov, drhow
673      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  rho_i, rho_j, rho_k
674      REAL(wp), DIMENSION(jpi,jpj)   ::  zcpx, zcpy    !W/D pressure filter
675      !!----------------------------------------------------------------------
676      !
677      !
678      !
679      IF( ln_wd ) THEN
680        DO jj = 2, jpjm1
681           DO ji = 2, jpim1 
682             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
683                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
684                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
685                  &                                                         > rn_wdmin1 + rn_wdmin2
686             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
687                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
688                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
689
690             IF(ll_tmp1) THEN
691               zcpx(ji,jj) = 1.0_wp
692             ELSE IF(ll_tmp2) THEN
693               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
694               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
695                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
696             ELSE
697               zcpx(ji,jj) = 0._wp
698             END IF
699     
700             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
701                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
702                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
703                  &                                                         > rn_wdmin1 + rn_wdmin2
704             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
705                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
706                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
707
708             IF(ll_tmp1) THEN
709               zcpy(ji,jj) = 1.0_wp
710             ELSE IF(ll_tmp2) THEN
711               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
712               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
713                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
714             ELSE
715               zcpy(ji,jj) = 0._wp
716             END IF
717           END DO
718        END DO
719        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
720      END IF
721
722      IF( kt == nit000 ) THEN
723         IF(lwp) WRITE(numout,*)
724         IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
725         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, density Jacobian with cubic polynomial scheme'
726      ENDIF
727
728      ! Local constant initialization
729      zcoef0 = - grav * 0.5_wp
730      z1_10  = 1._wp / 10._wp
731      z1_12  = 1._wp / 12._wp
732
733      !----------------------------------------------------------------------------------------
734      !  compute and store in provisional arrays elementary vertical and horizontal differences
735      !----------------------------------------------------------------------------------------
736
737!!bug gm   Not a true bug, but... dzz=e3w  for dzx, dzy verify what it is really
738
739      DO jk = 2, jpkm1
740         DO jj = 2, jpjm1
741            DO ji = fs_2, fs_jpim1   ! vector opt.
742               drhoz(ji,jj,jk) = rhd    (ji  ,jj  ,jk) - rhd    (ji,jj,jk-1)
743               dzz  (ji,jj,jk) = gde3w_n(ji  ,jj  ,jk) - gde3w_n(ji,jj,jk-1)
744               drhox(ji,jj,jk) = rhd    (ji+1,jj  ,jk) - rhd    (ji,jj,jk  )
745               dzx  (ji,jj,jk) = gde3w_n(ji+1,jj  ,jk) - gde3w_n(ji,jj,jk  )
746               drhoy(ji,jj,jk) = rhd    (ji  ,jj+1,jk) - rhd    (ji,jj,jk  )
747               dzy  (ji,jj,jk) = gde3w_n(ji  ,jj+1,jk) - gde3w_n(ji,jj,jk  )
748            END DO
749         END DO
750      END DO
751
752      !-------------------------------------------------------------------------
753      ! compute harmonic averages using eq. 5.18
754      !-------------------------------------------------------------------------
755      zep = 1.e-15
756
757!!bug  gm  drhoz not defined at level 1 and used (jk-1 with jk=2)
758!!bug  gm  idem for drhox, drhoy et ji=jpi and jj=jpj
759
760      DO jk = 2, jpkm1
761         DO jj = 2, jpjm1
762            DO ji = fs_2, fs_jpim1   ! vector opt.
763               cffw = 2._wp * drhoz(ji  ,jj  ,jk) * drhoz(ji,jj,jk-1)
764
765               cffu = 2._wp * drhox(ji+1,jj  ,jk) * drhox(ji,jj,jk  )
766               cffx = 2._wp * dzx  (ji+1,jj  ,jk) * dzx  (ji,jj,jk  )
767
768               cffv = 2._wp * drhoy(ji  ,jj+1,jk) * drhoy(ji,jj,jk  )
769               cffy = 2._wp * dzy  (ji  ,jj+1,jk) * dzy  (ji,jj,jk  )
770
771               IF( cffw > zep) THEN
772                  drhow(ji,jj,jk) = 2._wp *   drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1)   &
773                     &                    / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
774               ELSE
775                  drhow(ji,jj,jk) = 0._wp
776               ENDIF
777
778               dzw(ji,jj,jk) = 2._wp *   dzz(ji,jj,jk) * dzz(ji,jj,jk-1)   &
779                  &                  / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
780
781               IF( cffu > zep ) THEN
782                  drhou(ji,jj,jk) = 2._wp *   drhox(ji+1,jj,jk) * drhox(ji,jj,jk)   &
783                     &                    / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
784               ELSE
785                  drhou(ji,jj,jk ) = 0._wp
786               ENDIF
787
788               IF( cffx > zep ) THEN
789                  dzu(ji,jj,jk) = 2._wp *   dzx(ji+1,jj,jk) * dzx(ji,jj,jk)   &
790                     &                  / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
791               ELSE
792                  dzu(ji,jj,jk) = 0._wp
793               ENDIF
794
795               IF( cffv > zep ) THEN
796                  drhov(ji,jj,jk) = 2._wp *   drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk)   &
797                     &                    / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
798               ELSE
799                  drhov(ji,jj,jk) = 0._wp
800               ENDIF
801
802               IF( cffy > zep ) THEN
803                  dzv(ji,jj,jk) = 2._wp *   dzy(ji,jj+1,jk) * dzy(ji,jj,jk)   &
804                     &                  / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
805               ELSE
806                  dzv(ji,jj,jk) = 0._wp
807               ENDIF
808
809            END DO
810         END DO
811      END DO
812
813      !----------------------------------------------------------------------------------
814      ! apply boundary conditions at top and bottom using 5.36-5.37
815      !----------------------------------------------------------------------------------
816      drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:,  1  ) ) - 0.5_wp * drhow(:,:,  2  )
817      drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:,  1  ) ) - 0.5_wp * drhou(:,:,  2  )
818      drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:,  1  ) ) - 0.5_wp * drhov(:,:,  2  )
819
820      drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
821      drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
822      drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
823
824
825      !--------------------------------------------------------------
826      ! Upper half of top-most grid box, compute and store
827      !-------------------------------------------------------------
828
829!!bug gm   :  e3w-gde3w = 0.5*e3w  ....  and gde3w(2)-gde3w(1)=e3w(2) ....   to be verified
830!          true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
831
832      DO jj = 2, jpjm1
833         DO ji = fs_2, fs_jpim1   ! vector opt.
834            rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) )               &
835               &                   * (  rhd(ji,jj,1)                                     &
836               &                     + 0.5_wp * ( rhd    (ji,jj,2) - rhd    (ji,jj,1) )  &
837               &                              * ( e3w_n  (ji,jj,1) - gde3w_n(ji,jj,1) )  &
838               &                              / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) )  )
839         END DO
840      END DO
841
842!!bug gm    : here also, simplification is possible
843!!bug gm    : optimisation: 1/10 and 1/12 the division should be done before the loop
844
845      DO jk = 2, jpkm1
846         DO jj = 2, jpjm1
847            DO ji = fs_2, fs_jpim1   ! vector opt.
848
849               rho_k(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj,jk) + rhd    (ji,jj,jk-1) )                                   &
850                  &                     * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) )                                   &
851                  &            - grav * z1_10 * (                                                                   &
852                  &     ( drhow  (ji,jj,jk) - drhow  (ji,jj,jk-1) )                                                     &
853                  &   * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw  (ji,jj,jk) + dzw  (ji,jj,jk-1) ) )   &
854                  &   - ( dzw    (ji,jj,jk) - dzw    (ji,jj,jk-1) )                                                     &
855                  &   * ( rhd    (ji,jj,jk) - rhd    (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) )   &
856                  &                             )
857
858               rho_i(ji,jj,jk) = zcoef0 * ( rhd    (ji+1,jj,jk) + rhd    (ji,jj,jk) )                                   &
859                  &                     * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) )                                   &
860                  &            - grav* z1_10 * (                                                                    &
861                  &     ( drhou  (ji+1,jj,jk) - drhou  (ji,jj,jk) )                                                     &
862                  &   * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu  (ji+1,jj,jk) + dzu  (ji,jj,jk) ) )   &
863                  &   - ( dzu    (ji+1,jj,jk) - dzu    (ji,jj,jk) )                                                     &
864                  &   * ( rhd    (ji+1,jj,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) )   &
865                  &                            )
866
867               rho_j(ji,jj,jk) = zcoef0 * ( rhd    (ji,jj+1,jk) + rhd    (ji,jj,jk) )                                 &
868                  &                     * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) )                                   &
869                  &            - grav* z1_10 * (                                                                    &
870                  &     ( drhov  (ji,jj+1,jk) - drhov  (ji,jj,jk) )                                                     &
871                  &   * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv  (ji,jj+1,jk) + dzv  (ji,jj,jk) ) )   &
872                  &   - ( dzv    (ji,jj+1,jk) - dzv    (ji,jj,jk) )                                                     &
873                  &   * ( rhd    (ji,jj+1,jk) - rhd    (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) )   &
874                  &                            )
875
876            END DO
877         END DO
878      END DO
879      CALL lbc_lnk(rho_k,'W',1.)
880      CALL lbc_lnk(rho_i,'U',1.)
881      CALL lbc_lnk(rho_j,'V',1.)
882
883
884      ! ---------------
885      !  Surface value
886      ! ---------------
887      DO jj = 2, jpjm1
888         DO ji = fs_2, fs_jpim1   ! vector opt.
889            zhpi(ji,jj,1) = ( rho_k(ji+1,jj  ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj)
890            zhpj(ji,jj,1) = ( rho_k(ji  ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj)
891            IF( ln_wd ) THEN
892              zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
893              zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 
894            ENDIF
895            ! add to the general momentum trend
896            ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
897            va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
898         END DO
899      END DO
900
901      ! ----------------
902      !  interior value   (2=<jk=<jpkm1)
903      ! ----------------
904      DO jk = 2, jpkm1
905         DO jj = 2, jpjm1
906            DO ji = fs_2, fs_jpim1   ! vector opt.
907               ! hydrostatic pressure gradient along s-surfaces
908               zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1)                                &
909                  &           + (  ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk  ) )    &
910                  &              - ( rho_i(ji  ,jj,jk) - rho_i(ji,jj,jk-1) )  ) * r1_e1u(ji,jj)
911               zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1)                                &
912                  &           + (  ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk  ) )    &
913                  &               -( rho_j(ji,jj  ,jk) - rho_j(ji,jj,jk-1) )  ) * r1_e2v(ji,jj)
914               IF( ln_wd ) THEN
915                 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
916                 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 
917               ENDIF
918               ! add to the general momentum trend
919               ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
920               va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
921            END DO
922         END DO
923      END DO
924      !
925      !
926   END SUBROUTINE hpg_djc
927
928
929   SUBROUTINE hpg_prj( kt )
930      !!---------------------------------------------------------------------
931      !!                  ***  ROUTINE hpg_prj  ***
932      !!
933      !! ** Method  :   s-coordinate case.
934      !!      A Pressure-Jacobian horizontal pressure gradient method
935      !!      based on the constrained cubic-spline interpolation for
936      !!      all vertical coordinate systems
937      !!
938      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
939      !!----------------------------------------------------------------------
940      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear
941      INTEGER, INTENT(in) ::   kt                   ! ocean time-step index
942      !!
943      INTEGER  ::   ji, jj, jk, jkk                 ! dummy loop indices
944      REAL(wp) ::   zcoef0, znad                    ! local scalars
945      !
946      !! The local variables for the correction term
947      INTEGER  :: jk1, jis, jid, jjs, jjd
948      LOGICAL  :: ll_tmp1, ll_tmp2                  ! local logical variables
949      REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
950      REAL(wp) :: zrhdt1
951      REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
952      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zdept, zrhh
953      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
954      REAL(wp), DIMENSION(jpi,jpj)   ::   zsshu_n, zsshv_n
955      REAL(wp), DIMENSION(jpi,jpj)   ::  zcpx, zcpy    !W/D pressure filter
956      !!----------------------------------------------------------------------
957      !
958      !
959      IF( kt == nit000 ) THEN
960         IF(lwp) WRITE(numout,*)
961         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
962         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian'
963      ENDIF
964
965      ! Local constant initialization
966      zcoef0 = - grav
967      znad = 1._wp
968      IF( ln_linssh )   znad = 0._wp
969
970      IF( ln_wd ) THEN
971        DO jj = 2, jpjm1
972           DO ji = 2, jpim1 
973             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
974                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) .AND.            &
975                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji+1,jj) + ht_wd(ji+1,jj) ) &
976                  &                                                         > rn_wdmin1 + rn_wdmin2
977             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji+1,jj) ) > 1.E-12 ) .AND. (             &
978                  &    MAX(   sshn(ji,jj)               ,   sshn(ji+1,jj) ) >                &
979                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
980
981             IF(ll_tmp1) THEN
982               zcpx(ji,jj) = 1.0_wp
983             ELSE IF(ll_tmp2) THEN
984               ! no worries about  sshn(ji+1,jj) -  sshn(ji  ,jj) = 0, it won't happen ! here
985               zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_wd(ji+1,jj) - sshn(ji,jj) - ht_wd(ji,jj)) &
986                           &    / (sshn(ji+1,jj) -  sshn(ji  ,jj)) )
987             ELSE
988               zcpx(ji,jj) = 0._wp
989             END IF
990     
991             ll_tmp1 = MIN(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
992                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) .AND.            &
993                  &    MAX(   sshn(ji,jj) + ht_wd(ji,jj),   sshn(ji,jj+1) + ht_wd(ji,jj+1) ) &
994                  &                                                         > rn_wdmin1 + rn_wdmin2
995             ll_tmp2 = ( ABS( sshn(ji,jj)               -   sshn(ji,jj+1) ) > 1.E-12 ) .AND. (             &
996                  &    MAX(   sshn(ji,jj)               ,   sshn(ji,jj+1) ) >                &
997                  &    MAX( -ht_wd(ji,jj)               , -ht_wd(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
998
999             IF(ll_tmp1) THEN
1000               zcpy(ji,jj) = 1.0_wp
1001             ELSE IF(ll_tmp2) THEN
1002               ! no worries about  sshn(ji,jj+1) -  sshn(ji,jj  ) = 0, it won't happen ! here
1003               zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_wd(ji,jj+1) - sshn(ji,jj) - ht_wd(ji,jj)) &
1004                           &    / (sshn(ji,jj+1) -  sshn(ji,jj  )) )
1005             ELSE
1006               zcpy(ji,jj) = 0._wp
1007             END IF
1008           END DO
1009        END DO
1010        CALL lbc_lnk( zcpx, 'U', 1._wp )    ;   CALL lbc_lnk( zcpy, 'V', 1._wp )
1011      END IF
1012
1013      ! Clean 3-D work arrays
1014      zhpi(:,:,:) = 0._wp
1015      zrhh(:,:,:) = rhd(:,:,:)
1016
1017      ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
1018      DO jj = 1, jpj
1019        DO ji = 1, jpi
1020          jk = mbkt(ji,jj)+1
1021          IF(     jk <=  0   ) THEN   ;   zrhh(ji,jj,    :   ) = 0._wp
1022          ELSEIF( jk ==  1   ) THEN   ;   zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)
1023          ELSEIF( jk < jpkm1 ) THEN
1024             DO jkk = jk+1, jpk
1025                zrhh(ji,jj,jkk) = interp1(gde3w_n(ji,jj,jkk  ), gde3w_n(ji,jj,jkk-1),   &
1026                   &                      gde3w_n(ji,jj,jkk-2), rhd    (ji,jj,jkk-1), rhd(ji,jj,jkk-2))
1027             END DO
1028          ENDIF
1029        END DO
1030      END DO
1031
1032      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
1033      DO jj = 1, jpj
1034         DO ji = 1, jpi
1035            zdept(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) - sshn(ji,jj) * znad
1036         END DO
1037      END DO
1038
1039      DO jk = 2, jpk
1040         DO jj = 1, jpj
1041            DO ji = 1, jpi
1042               zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w_n(ji,jj,jk)
1043            END DO
1044         END DO
1045      END DO
1046
1047      fsp(:,:,:) = zrhh (:,:,:)
1048      xsp(:,:,:) = zdept(:,:,:)
1049
1050      ! Construct the vertical density profile with the
1051      ! constrained cubic spline interpolation
1052      ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3
1053      CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
1054
1055      ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
1056      DO jj = 2, jpj
1057        DO ji = 2, jpi
1058          zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1),  &
1059             &                                              csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1)
1060
1061          ! assuming linear profile across the top half surface layer
1062          zhpi(ji,jj,1) =  0.5_wp * e3w_n(ji,jj,1) * zrhdt1
1063        END DO
1064      END DO
1065
1066      ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
1067      DO jk = 2, jpkm1
1068        DO jj = 2, jpj
1069          DO ji = 2, jpi
1070            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) +                                  &
1071               &             integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk),   &
1072               &                           asp  (ji,jj,jk-1), bsp  (ji,jj,jk-1), &
1073               &                           csp  (ji,jj,jk-1), dsp  (ji,jj,jk-1)  )
1074          END DO
1075        END DO
1076      END DO
1077
1078      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
1079
1080      ! Prepare zsshu_n and zsshv_n
1081      DO jj = 2, jpjm1
1082        DO ji = 2, jpim1
1083!!gm BUG ?    if it is ssh at u- & v-point then it should be:
1084!          zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * &
1085!                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
1086!          zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * &
1087!                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
1088!!gm not this:
1089          zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &
1090                         & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 
1091          zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &
1092                         & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 
1093        END DO
1094      END DO
1095
1096      CALL lbc_lnk (zsshu_n, 'U', 1.)
1097      CALL lbc_lnk (zsshv_n, 'V', 1.)
1098
1099      DO jj = 2, jpjm1
1100        DO ji = 2, jpim1
1101          zu(ji,jj,1) = - ( e3u_n(ji,jj,1) - zsshu_n(ji,jj) * znad) 
1102          zv(ji,jj,1) = - ( e3v_n(ji,jj,1) - zsshv_n(ji,jj) * znad)
1103        END DO
1104      END DO
1105
1106      DO jk = 2, jpkm1
1107        DO jj = 2, jpjm1
1108          DO ji = 2, jpim1
1109            zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk)
1110            zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk)
1111          END DO
1112        END DO
1113      END DO
1114
1115      DO jk = 1, jpkm1
1116        DO jj = 2, jpjm1
1117          DO ji = 2, jpim1
1118            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u_n(ji,jj,jk)
1119            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v_n(ji,jj,jk)
1120          END DO
1121        END DO
1122      END DO
1123
1124      DO jk = 1, jpkm1
1125        DO jj = 2, jpjm1
1126          DO ji = 2, jpim1
1127            zu(ji,jj,jk) = MIN(  zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
1128            zu(ji,jj,jk) = MAX(  zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) )  )
1129            zv(ji,jj,jk) = MIN(  zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
1130            zv(ji,jj,jk) = MAX(  zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) )  )
1131          END DO
1132        END DO
1133      END DO
1134
1135
1136      DO jk = 1, jpkm1
1137        DO jj = 2, jpjm1
1138          DO ji = 2, jpim1
1139            zpwes = 0._wp; zpwed = 0._wp
1140            zpnss = 0._wp; zpnsd = 0._wp
1141            zuijk = zu(ji,jj,jk)
1142            zvijk = zv(ji,jj,jk)
1143
1144            !!!!!     for u equation
1145            IF( jk <= mbku(ji,jj) ) THEN
1146               IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
1147                 jis = ji + 1; jid = ji
1148               ELSE
1149                 jis = ji;     jid = ji +1
1150               ENDIF
1151
1152               ! integrate the pressure on the shallow side
1153               jk1 = jk
1154               DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
1155                 IF( jk1 == mbku(ji,jj) ) THEN
1156                   zuijk = -zdept(jis,jj,jk1)
1157                   EXIT
1158                 ENDIF
1159                 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
1160                 zpwes = zpwes +                                    &
1161                      integ_spline(zdept(jis,jj,jk1), zdeps,            &
1162                             asp(jis,jj,jk1),    bsp(jis,jj,jk1), &
1163                             csp(jis,jj,jk1),    dsp(jis,jj,jk1))
1164                 jk1 = jk1 + 1
1165               END DO
1166
1167               ! integrate the pressure on the deep side
1168               jk1 = jk
1169               DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
1170                 IF( jk1 == 1 ) THEN
1171                   zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
1172                   zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
1173                                                     bsp(jid,jj,1),   csp(jid,jj,1), &
1174                                                     dsp(jid,jj,1)) * zdeps
1175                   zpwed  = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
1176                   EXIT
1177                 ENDIF
1178                 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
1179                 zpwed = zpwed +                                        &
1180                        integ_spline(zdeps,              zdept(jid,jj,jk1), &
1181                               asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1),  &
1182                               csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
1183                 jk1 = jk1 - 1
1184               END DO
1185
1186               ! update the momentum trends in u direction
1187
1188               zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) )
1189               IF( .NOT.ln_linssh ) THEN
1190                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
1191                    &    ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )
1192                ELSE
1193                 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
1194               ENDIF
1195               IF( ln_wd ) THEN
1196                  zdpdx1 = zdpdx1 * zcpx(ji,jj)
1197                  zdpdx2 = zdpdx2 * zcpx(ji,jj)
1198               ENDIF
1199               ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 
1200            ENDIF
1201
1202            !!!!!     for v equation
1203            IF( jk <= mbkv(ji,jj) ) THEN
1204               IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
1205                 jjs = jj + 1; jjd = jj
1206               ELSE
1207                 jjs = jj    ; jjd = jj + 1
1208               ENDIF
1209
1210               ! integrate the pressure on the shallow side
1211               jk1 = jk
1212               DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
1213                 IF( jk1 == mbkv(ji,jj) ) THEN
1214                   zvijk = -zdept(ji,jjs,jk1)
1215                   EXIT
1216                 ENDIF
1217                 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
1218                 zpnss = zpnss +                                      &
1219                        integ_spline(zdept(ji,jjs,jk1), zdeps,            &
1220                               asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), &
1221                               csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) )
1222                 jk1 = jk1 + 1
1223               END DO
1224
1225               ! integrate the pressure on the deep side
1226               jk1 = jk
1227               DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
1228                 IF( jk1 == 1 ) THEN
1229                   zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
1230                   zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
1231                                                     bsp(ji,jjd,1),   csp(ji,jjd,1), &
1232                                                     dsp(ji,jjd,1) ) * zdeps
1233                   zpnsd  = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
1234                   EXIT
1235                 ENDIF
1236                 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
1237                 zpnsd = zpnsd +                                        &
1238                        integ_spline(zdeps,              zdept(ji,jjd,jk1), &
1239                               asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
1240                               csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
1241                 jk1 = jk1 - 1
1242               END DO
1243
1244
1245               ! update the momentum trends in v direction
1246
1247               zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )
1248               IF( .NOT.ln_linssh ) THEN
1249                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
1250                          ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )
1251               ELSE
1252                  zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
1253               ENDIF
1254               IF( ln_wd ) THEN
1255                  zdpdy1 = zdpdy1 * zcpy(ji,jj)
1256                  zdpdy2 = zdpdy2 * zcpy(ji,jj)
1257               ENDIF
1258
1259               va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
1260            ENDIF
1261               !
1262            END DO
1263         END DO
1264      END DO
1265      !
1266      !
1267   END SUBROUTINE hpg_prj
1268
1269
1270   SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
1271      !!----------------------------------------------------------------------
1272      !!                 ***  ROUTINE cspline  ***
1273      !!
1274      !! ** Purpose :   constrained cubic spline interpolation
1275      !!
1276      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3
1277      !!
1278      !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
1279      !!----------------------------------------------------------------------
1280      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   fsp, xsp           ! value and coordinate
1281      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   asp, bsp, csp, dsp ! coefficients of the interpoated function
1282      INTEGER                   , INTENT(in   ) ::   polynomial_type    ! 1: cubic spline   ;   2: Linear
1283      !
1284      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
1285      INTEGER  ::   jpi, jpj, jpkm1
1286      REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
1287      REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha
1288      REAL(wp) ::   zdf(size(fsp,3))
1289      !!----------------------------------------------------------------------
1290      !
1291!!gm  WHAT !!!!!   THIS IS VERY DANGEROUS !!!!! 
1292      jpi   = size(fsp,1)
1293      jpj   = size(fsp,2)
1294      jpkm1 = MAX( 1, size(fsp,3) - 1 )
1295      !
1296      IF (polynomial_type == 1) THEN     ! Constrained Cubic Spline
1297         DO ji = 1, jpi
1298            DO jj = 1, jpj
1299           !!Fritsch&Butland's method, 1984 (preferred, but more computation)
1300           !    DO jk = 2, jpkm1-1
1301           !       zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)
1302           !       zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1303           !       zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1
1304           !       zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2
1305           !
1306           !       zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
1307           !
1308           !       IF(zdf1 * zdf2 <= 0._wp) THEN
1309           !           zdf(jk) = 0._wp
1310           !       ELSE
1311           !         zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
1312           !       ENDIF
1313           !    END DO
1314
1315           !!Simply geometric average
1316               DO jk = 2, jpkm1-1
1317                  zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1))
1318                  zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk))
1319
1320                  IF(zdf1 * zdf2 <= 0._wp) THEN
1321                     zdf(jk) = 0._wp
1322                  ELSE
1323                     zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
1324                  ENDIF
1325               END DO
1326
1327               zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
1328                          &          ( xsp(ji,jj,2) - xsp(ji,jj,1) )           -  0.5_wp * zdf(2)
1329               zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
1330                          &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1)
1331
1332               DO jk = 1, jpkm1 - 1
1333                 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1334                 ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
1335                 ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
1336                 zddf1  = -2._wp * ztmp1 + ztmp2
1337                 ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
1338                 zddf2  =  2._wp * ztmp1 - ztmp2
1339
1340                 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
1341                 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
1342                 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
1343                               & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
1344                               & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
1345                               &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
1346                 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
1347                               &                (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
1348                               &                 dsp(ji,jj,jk) * xsp(ji,jj,jk))))
1349               END DO
1350            END DO
1351         END DO
1352
1353      ELSEIF ( polynomial_type == 2 ) THEN     ! Linear
1354         DO ji = 1, jpi
1355            DO jj = 1, jpj
1356               DO jk = 1, jpkm1-1
1357                  zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1358                  ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
1359
1360                  dsp(ji,jj,jk) = 0._wp
1361                  csp(ji,jj,jk) = 0._wp
1362                  bsp(ji,jj,jk) = ztmp1 / zdxtmp
1363                  asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
1364               END DO
1365            END DO
1366         END DO
1367
1368      ELSE
1369           CALL ctl_stop( 'invalid polynomial type in cspline' )
1370      ENDIF
1371
1372   END SUBROUTINE cspline
1373
1374
1375   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f)
1376      !!----------------------------------------------------------------------
1377      !!                 ***  ROUTINE interp1  ***
1378      !!
1379      !! ** Purpose :   1-d linear interpolation
1380      !!
1381      !! ** Method  :   interpolation is straight forward
1382      !!                extrapolation is also permitted (no value limit)
1383      !!----------------------------------------------------------------------
1384      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr
1385      REAL(wp)             ::  f ! result of the interpolation (extrapolation)
1386      REAL(wp)             ::  zdeltx
1387      !!----------------------------------------------------------------------
1388      !
1389      zdeltx = xr - xl
1390      IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN
1391         f = 0.5_wp * (fl + fr)
1392      ELSE
1393         f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
1394      ENDIF
1395      !
1396   END FUNCTION interp1
1397
1398
1399   FUNCTION interp2( x, a, b, c, d )  RESULT(f)
1400      !!----------------------------------------------------------------------
1401      !!                 ***  ROUTINE interp1  ***
1402      !!
1403      !! ** Purpose :   1-d constrained cubic spline interpolation
1404      !!
1405      !! ** Method  :  cubic spline interpolation
1406      !!
1407      !!----------------------------------------------------------------------
1408      REAL(wp), INTENT(in) ::  x, a, b, c, d
1409      REAL(wp)             ::  f ! value from the interpolation
1410      !!----------------------------------------------------------------------
1411      !
1412      f = a + x* ( b + x * ( c + d * x ) )
1413      !
1414   END FUNCTION interp2
1415
1416
1417   FUNCTION interp3( x, a, b, c, d )  RESULT(f)
1418      !!----------------------------------------------------------------------
1419      !!                 ***  ROUTINE interp1  ***
1420      !!
1421      !! ** Purpose :   Calculate the first order of deriavtive of
1422      !!                a cubic spline function y=a+b*x+c*x^2+d*x^3
1423      !!
1424      !! ** Method  :   f=dy/dx=b+2*c*x+3*d*x^2
1425      !!
1426      !!----------------------------------------------------------------------
1427      REAL(wp), INTENT(in) ::  x, a, b, c, d
1428      REAL(wp)             ::  f ! value from the interpolation
1429      !!----------------------------------------------------------------------
1430      !
1431      f = b + x * ( 2._wp * c + 3._wp * d * x)
1432      !
1433   END FUNCTION interp3
1434
1435
1436   FUNCTION integ_spline( xl, xr, a, b, c, d )  RESULT(f)
1437      !!----------------------------------------------------------------------
1438      !!                 ***  ROUTINE interp1  ***
1439      !!
1440      !! ** Purpose :   1-d constrained cubic spline integration
1441      !!
1442      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr
1443      !!
1444      !!----------------------------------------------------------------------
1445      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d
1446      REAL(wp)             ::  za1, za2, za3
1447      REAL(wp)             ::  f                   ! integration result
1448      !!----------------------------------------------------------------------
1449      !
1450      za1 = 0.5_wp * b
1451      za2 = c / 3.0_wp
1452      za3 = 0.25_wp * d
1453      !
1454      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
1455         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
1456      !
1457   END FUNCTION integ_spline
1458
1459   !!======================================================================
1460END MODULE dynhpg
1461
Note: See TracBrowser for help on using the repository browser.