source: branches/UKMO/dev_3.6_FVPS/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 8476

Last change on this file since 8476 was 8476, checked in by davestorkey, 3 years ago

UKMO/dev_3.6_FVPS : implement FVPS HPG scheme.

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