source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 18 months ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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