New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
scoord_gen.F90 in branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src – NEMO

source: branches/2015/dev_r5187_UKMO13_simplification/NEMOGCM/TOOLS/SCOORD_GEN/src/scoord_gen.F90 @ 5255

Last change on this file since 5255 was 5255, checked in by timgraham, 9 years ago

First attempt at making a tool to move Scoord generation offline.
File size: 39.9 KB
Line 
1PROGRAM SCOORD_GEN
2      !!----------------------------------------------------------------------
3      !!                  ***  PROGRAM scoord_gen  ***
4      !!                     
5      !! ** Purpose :   define the s-coordinate system
6      !!
7      !! ** Method  :   s-coordinate
8      !!         The depth of model levels is defined as the product of an
9      !!      analytical function by the local bathymetry, while the vertical
10      !!      scale factors are defined as the product of the first derivative
11      !!      of the analytical function by the bathymetry.
12      !!      (this solution save memory as depth and scale factors are not
13      !!      3d fields)
14      !!          - Read bathymetry (in meters) at t-point and compute the
15      !!         bathymetry at u-, v-, and f-points.
16      !!            hbatu = mi( hbatt )
17      !!            hbatv = mj( hbatt )
18      !!            hbatf = mi( mj( hbatt ) )
19      !!          - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
20      !!         function and its derivative given as function.
21      !!            z_gsigt(k) = fssig (k    )
22      !!            z_gsigw(k) = fssig (k-0.5)
23      !!            z_esigt(k) = fsdsig(k    )
24      !!            z_esigw(k) = fsdsig(k-0.5)
25      !!      Three options for stretching are give, and they can be modified
26      !!      following the users requirements. Nevertheless, the output as
27      !!      well as the way to compute the model levels and scale factors
28      !!      must be respected in order to insure second order accuracy
29      !!      schemes.
30      !!
31      !!      The three methods for stretching available are:
32      !!
33      !!           s_sh94 (Song and Haidvogel 1994)
34      !!                a sinh/tanh function that allows sigma and stretched sigma
35      !!
36      !!           s_sf12 (Siddorn and Furner 2012?)
37      !!                allows the maintenance of fixed surface and or
38      !!                bottom cell resolutions (cf. geopotential coordinates)
39      !!                within an analytically derived stretched S-coordinate framework.
40      !!
41      !!          s_tanh  (Madec et al 1996)
42      !!                a cosh/tanh function that gives stretched coordinates       
43      !!
44      !!----------------------------------------------------------------------
45      !
46      IMPLICIT NONE
47
48      USE utils
49
50     !!----------------------------------------------------------------------
51      !
52      !
53      OPEN( UNIT=numnam, FILE='namelist', FORM='FORMATTED', STATUS='OLD' )
54      READ( numnam, namzgr_sco )
55      CLOSE( numnam )
56      jpk = rn_jpk
57
58      WRITE(numout,*)
59      WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
60      WRITE(numout,*) '~~~~~~~~~~~'
61      WRITE(numout,*) '   Namelist namzgr_sco'
62      WRITE(numout,*) '     stretching coeffs '
63      WRITE(numout,*) '        Number of levels                             rn_jpk        = ',rn_jpk
64      WRITE(numout,*) '        maximum depth of s-bottom surface (>0)       rn_sbot_max   = ',rn_sbot_max
65      WRITE(numout,*) '        minimum depth of s-bottom surface (>0)       rn_sbot_min   = ',rn_sbot_min
66      WRITE(numout,*) '        Critical depth                               rn_hc         = ',rn_hc
67      WRITE(numout,*) '        maximum cut-off r-value allowed              rn_rmax       = ',rn_rmax
68      WRITE(numout,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
69      WRITE(numout,*) '        Song and Haidvogel 1994 stretching coefficients'
70      WRITE(numout,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
71      WRITE(numout,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
72      WRITE(numout,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
73      WRITE(numout,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
74      WRITE(numout,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
75      WRITE(numout,*) '        Siddorn and Furner 2012 stretching coefficients'
76      WRITE(numout,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
77      WRITE(numout,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
78      WRITE(numout,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
79      WRITE(numout,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
80      WRITE(numout,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
81      WRITE(numout,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
82
83     
84      ! Read in bathy, jpi and jpj from bathy.nc
85      CALL read_bathy()
86
87      !Allocate all other allocatable variables
88      ios = dom_oce_alloc()
89      IF( ios .NE. 0) THEN
90         WRITE(0,*) 'Unable to allocate all arrays'
91         STOP 1
92      ENDIF
93
94      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
95      hifu(:,:) = rn_sbot_min
96      hifv(:,:) = rn_sbot_min
97      hiff(:,:) = rn_sbot_min
98
99      !                                        ! set maximum ocean depth
100      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
101
102      DO jj = 1, jpj
103         DO ji = 1, jpi
104           IF( bathy(ji,jj) > 0._wp )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
105         END DO
106      END DO
107      !                                        ! =============================
108      !                                        ! Define the envelop bathymetry   (hbatt)
109      !                                        ! =============================
110      ! use r-value to create hybrid coordinates
111      zenv(:,:) = bathy(:,:)
112      !
113      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
114      DO jj = 1, jpj
115         DO ji = 1, jpi
116           IF( bathy(ji,jj) == 0._wp ) THEN
117             iip1 = MIN( ji+1, jpi )
118             ijp1 = MIN( jj+1, jpj )
119             iim1 = MAX( ji-1, 1 )
120             ijm1 = MAX( jj-1, 1 )
121             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
122        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
123               zenv(ji,jj) = rn_sbot_min
124             ENDIF
125           ENDIF
126         END DO
127      END DO
128      !
129      ! smooth the bathymetry (if required)
130      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea)
131      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
132      !
133      jl = 0
134      zrmax = 1._wp
135      !   
136      !     
137      ! set scaling factor used in reducing vertical gradients
138      zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
139      !
140      ! initialise temporary evelope depth arrays
141      ztmpi1(:,:) = zenv(:,:)
142      ztmpi2(:,:) = zenv(:,:)
143      ztmpj1(:,:) = zenv(:,:)
144      ztmpj2(:,:) = zenv(:,:)
145      !
146      ! initialise temporary r-value arrays
147      zri(:,:) = 1._wp
148      zrj(:,:) = 1._wp
149      !                                                            ! ================ !
150      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) !  Iterative loop  !
151         !                                                         ! ================ !
152         jl = jl + 1
153         zrmax = 0._wp
154         ! we set zrmax from previous r-values (zri and zrj) first
155         ! if set after current r-value calculation (as previously)
156         ! we could exit DO WHILE prematurely before checking r-value
157         ! of current zenv
158         DO jj = 1, nlcj
159            DO ji = 1, nlci
160               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
161            END DO
162         END DO
163         zri(:,:) = 0._wp
164         zrj(:,:) = 0._wp
165         DO jj = 1, nlcj
166            DO ji = 1, nlci
167               iip1 = MIN( ji+1, nlci )      ! force zri = 0 on last line (ji=ncli+1 to jpi)
168               ijp1 = MIN( jj+1, nlcj )      ! force zrj = 0 on last raw  (jj=nclj+1 to jpj)
169               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
170                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
171               END IF
172               IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
173                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
174               END IF
175               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
176               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
177               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
178               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
179            END DO
180         END DO
181         !
182         WRITE(*,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
183         !
184         DO jj = 1, nlcj
185            DO ji = 1, nlci
186               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
187            END DO
188         END DO
189         ! apply lateral boundary condition   CAUTION: keep the value when the lbc field is zero
190         !                                                  ! ================ !
191      END DO                                                !     End loop     !
192      !                                                     ! ================ !
193      DO jj = 1, jpj
194         DO ji = 1, jpi
195            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
196         END DO
197      END DO
198      !
199      ! Envelope bathymetry saved in hbatt
200      hbatt(:,:) = zenv(:,:) 
201      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
202         CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
203         DO jj = 1, jpj
204            DO ji = 1, jpi
205               ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
206               hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
207            END DO
208         END DO
209      ENDIF
210      !
211      !                                        ! ==============================
212      !                                        !   hbatu, hbatv, hbatf fields
213      !                                        ! ==============================
214      WRITE(numout,*)
215      WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
216      hbatu(:,:) = rn_sbot_min
217      hbatv(:,:) = rn_sbot_min
218      hbatf(:,:) = rn_sbot_min
219      DO jj = 1, jpjm1
220        DO ji = 1, jpim1   ! NO vector opt.
221           hbatu(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
222           hbatv(ji,jj) = 0.50_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
223           hbatf(ji,jj) = 0.25_wp * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
224              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
225        END DO
226      END DO
227      !
228      zhbat(:,:) = hbatu(:,:)   
229      DO jj = 1, jpj
230         DO ji = 1, jpi
231            IF( hbatu(ji,jj) == 0._wp ) THEN
232               IF( zhbat(ji,jj) == 0._wp )   hbatu(ji,jj) = rn_sbot_min
233               IF( zhbat(ji,jj) /= 0._wp )   hbatu(ji,jj) = zhbat(ji,jj)
234            ENDIF
235         END DO
236      END DO
237      zhbat(:,:) = hbatv(:,:) 
238      DO jj = 1, jpj
239         DO ji = 1, jpi
240            IF( hbatv(ji,jj) == 0._wp ) THEN
241               IF( zhbat(ji,jj) == 0._wp )   hbatv(ji,jj) = rn_sbot_min
242               IF( zhbat(ji,jj) /= 0._wp )   hbatv(ji,jj) = zhbat(ji,jj)
243            ENDIF
244         END DO
245      END DO
246      zhbat(:,:) = hbatf(:,:)   
247      DO jj = 1, jpj
248         DO ji = 1, jpi
249            IF( hbatf(ji,jj) == 0._wp ) THEN
250               IF( zhbat(ji,jj) == 0._wp )   hbatf(ji,jj) = rn_sbot_min
251               IF( zhbat(ji,jj) /= 0._wp )   hbatf(ji,jj) = zhbat(ji,jj)
252            ENDIF
253         END DO
254      END DO
255
256!!bug:  key_helsinki a verifer
257      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
258      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
259      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
260      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
261
262         WRITE(numout,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
263            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
264         WRITE(numout,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
265            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
266         WRITE(numout,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
267            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
268         WRITE(numout,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
269            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
270!! helsinki
271
272      !                                            ! =======================
273      !                                            !   s-coordinate fields     (gdep., e3.)
274      !                                            ! =======================
275      !
276      ! non-dimensional "sigma" for model level depth at w- and t-levels
277
278!========================================================================
279! Song and Haidvogel  1994 (ln_s_sh94=T)
280! Siddorn and Furner 2012 (ln_sf12=T)
281! or  tanh function       (both false)                   
282!========================================================================
283      IF      ( ln_s_sh94 ) THEN
284                           CALL s_sh94()
285      ELSE IF ( ln_s_sf12 ) THEN
286                           CALL s_sf12()
287      ELSE                 
288                           CALL s_tanh()
289      ENDIF 
290
291      !
292      where (e3t_0   (:,:,:).eq.0.0)  e3t_0(:,:,:) = 1.0
293      where (e3u_0   (:,:,:).eq.0.0)  e3u_0(:,:,:) = 1.0
294      where (e3v_0   (:,:,:).eq.0.0)  e3v_0(:,:,:) = 1.0
295      where (e3f_0   (:,:,:).eq.0.0)  e3f_0(:,:,:) = 1.0
296      where (e3w_0   (:,:,:).eq.0.0)  e3w_0(:,:,:) = 1.0
297      where (e3uw_0  (:,:,:).eq.0.0)  e3uw_0(:,:,:) = 1.0
298      where (e3vw_0  (:,:,:).eq.0.0)  e3vw_0(:,:,:) = 1.0
299
300#if defined key_agrif
301      ! Ensure meaningful vertical scale factors in ghost lines/columns
302      ! TODO: How can we deal with this in offline tool?
303      IF( .NOT. Agrif_Root() ) THEN
304         ! 
305         IF((nbondi == -1).OR.(nbondi == 2)) THEN
306            e3u_0(1,:,:) = e3u_0(2,:,:)
307         ENDIF
308         !
309         IF((nbondi ==  1).OR.(nbondi == 2)) THEN
310            e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:)
311         ENDIF
312         !
313         IF((nbondj == -1).OR.(nbondj == 2)) THEN
314            e3v_0(:,1,:) = e3v_0(:,2,:)
315         ENDIF
316         !
317         IF((nbondj ==  1).OR.(nbondj == 2)) THEN
318            e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:)
319         ENDIF
320         !
321      ENDIF
322#endif
323!!
324      ! HYBRID :
325      DO jj = 1, jpj
326         DO ji = 1, jpi
327            DO jk = 1, jpk-1
328               IF( scobot(ji,jj) >= fsdept(ji,jj,jk) )   mbathy(ji,jj) = MAX( 2, jk )
329            END DO
330            IF( scobot(ji,jj) == 0._wp               )   mbathy(ji,jj) = 0
331         END DO
332      END DO
333      WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
334
335      ! min max values over the domain
336      WRITE(numout,*) ' MIN val mbathy  ', MINVAL( mbathy(:,:)    ), ' MAX ', MAXVAL( mbathy(:,:) )
337      WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ),   &
338         &                          ' w ', MINVAL( gdepw_0(:,:,:) ), '3w '  , MINVAL( gdep3w_0(:,:,:) )
339      WRITE(numout,*) ' MIN val e3    t ', MINVAL( e3t_0  (:,:,:) ), ' f '  , MINVAL( e3f_0   (:,:,:) ),   &
340         &                          ' u ', MINVAL( e3u_0  (:,:,:) ), ' u '  , MINVAL( e3v_0   (:,:,:) ),   &
341         &                          ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw'  , MINVAL( e3vw_0  (:,:,:) ),   &
342         &                          ' w ', MINVAL( e3w_0  (:,:,:) )
343
344      WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ),   &
345         &                          ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w '  , MAXVAL( gdep3w_0(:,:,:) )
346      WRITE(numout,*) ' MAX val e3    t ', MAXVAL( e3t_0  (:,:,:) ), ' f '  , MAXVAL( e3f_0   (:,:,:) ),   &
347         &                          ' u ', MAXVAL( e3u_0  (:,:,:) ), ' u '  , MAXVAL( e3v_0   (:,:,:) ),   &
348         &                          ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw'  , MAXVAL( e3vw_0  (:,:,:) ),   &
349         &                          ' w ', MAXVAL( e3w_0  (:,:,:) )
350     
351         ! selected vertical profiles
352      WRITE(numout,*)
353      WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
354      WRITE(numout,*) ' ~~~~~~  --------------------'
355      WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
356      WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk),     &
357         &                                 e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
358      DO jj = mj0(20), mj1(20)
359         DO ji = mi0(20), mi1(20)
360            WRITE(numout,*)
361            WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
362            WRITE(numout,*) ' ~~~~~~  --------------------'
363            WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
364            WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
365               &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
366         END DO
367      END DO
368      DO jj = mj0(74), mj1(74)
369         DO ji = mi0(100), mi1(100)
370            WRITE(numout,*)
371            WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k)   bathy = ', bathy(ji,jj), hbatt(ji,jj)
372            WRITE(numout,*) ' ~~~~~~  --------------------'
373            WRITE(numout,"(9x,' level  gdept_0   gdepw_0   e3t_0    e3w_0')")
374            WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk),     &
375               &                                 e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
376         END DO
377      END DO
378
379!================================================================================
380! check the coordinate makes sense
381!================================================================================
382      DO ji = 1, jpi
383         DO jj = 1, jpj
384
385            IF( hbatt(ji,jj) > 0._wp) THEN
386               DO jk = 1, mbathy(ji,jj)
387                 ! check coordinate is monotonically increasing
388                 IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
389                    WRITE(numout,*) 'ERROR zgr_sco :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk
390                    WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
391                    WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
392                    CALL STOP 1
393                 ENDIF
394                 ! and check it has never gone negative
395                 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
396                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk
397                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
398                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
399                    CALL STOP 1
400                 ENDIF
401                 ! and check it never exceeds the total depth
402                 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
403                    WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk
404                    WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
405                    CALL STOP 1
406                 ENDIF
407               END DO
408
409               DO jk = 1, mbathy(ji,jj)-1
410                 ! and check it never exceeds the total depth
411                IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
412                    WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk
413                    WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
414                    CALL STOP 1
415                 ENDIF
416               END DO
417
418            ENDIF
419
420         END DO
421      END DO
422      !
423      !
424
425!!======================================================================
426   SUBROUTINE s_sh94()
427
428      !!----------------------------------------------------------------------
429      !!                  ***  ROUTINE s_sh94  ***
430      !!                     
431      !! ** Purpose :   stretch the s-coordinate system
432      !!
433      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
434      !!                mixed S/sigma coordinate
435      !!
436      !! Reference : Song and Haidvogel 1994.
437      !!----------------------------------------------------------------------
438      !
439      INTEGER  ::   ji, jj, jk           ! dummy loop argument
440      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
441      !
442      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
443      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 
444      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
445
446      ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) )
447      ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 
448      ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 
449      ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
450
451      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
452      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
453      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
454      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
455
456      DO ji = 1, jpi
457         DO jj = 1, jpj
458
459            IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
460               DO jk = 1, jpk
461                  z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
462                  z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp)       , rn_bb )
463               END DO
464            ELSE ! shallow water, uniform sigma
465               DO jk = 1, jpk
466                  z_gsigw3(ji,jj,jk) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
467                  z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
468                  END DO
469            ENDIF
470            !
471            DO jk = 1, jpk-1
472               z_esigt3(ji,jj,jk  ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
473               z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
474            END DO
475            z_esigw3(ji,jj,1  ) = 2._wp * ( z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ) )
476            z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
477            !
478            ! Coefficients for vertical depth as the sum of e3w scale factors
479            z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
480            DO jk = 2, jpk
481               z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
482            END DO
483            !
484            DO jk = 1, jpk
485               zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpk-1,wp)
486               zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp)
487               gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
488               gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
489               gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )
490            END DO
491           !
492         END DO   ! for all jj's
493      END DO    ! for all ji's
494
495      DO ji = 1, jpim1
496         DO jj = 1, jpjm1
497            DO jk = 1, jpk
498               z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) )   &
499                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
500               z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) )   &
501                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
502               z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk)     &
503                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) )   &
504                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
505               z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) )   &
506                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
507               z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) )   &
508                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
509               !
510               e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
511               e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
512               e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
513               e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
514               !
515               e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
516               e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
517               e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpk-1,wp) )
518            END DO
519        END DO
520      END DO
521
522      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3)
523      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
524
525   END SUBROUTINE s_sh94
526
527   SUBROUTINE s_sf12
528
529      !!----------------------------------------------------------------------
530      !!                  ***  ROUTINE s_sf12 ***
531      !!                     
532      !! ** Purpose :   stretch the s-coordinate system
533      !!
534      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
535      !!                mixed S/sigma/Z coordinate
536      !!
537      !!                This method allows the maintenance of fixed surface and or
538      !!                bottom cell resolutions (cf. geopotential coordinates)
539      !!                within an analytically derived stretched S-coordinate framework.
540      !!
541      !!
542      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
543      !!----------------------------------------------------------------------
544      !
545      INTEGER  ::   ji, jj, jk           ! dummy loop argument
546      REAL(wp) ::   zsmth               ! smoothing around critical depth
547      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
548      !
549      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
550      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3 
551      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
552      !
553      ALLOCATE( z_gsigw3(jpi,jpj,jpk), z_gsigt3(jpi,jpj,jpk), z_gsi3w3(jpi,jpj,jpk) )
554      ALLOCATE( z_esigt3(jpi,jpj,jpk), z_esigw3(jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk)) 
555      ALLOCATE( z_esigtv3(jpi,jpj,jpk), z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk)) 
556      ALLOCATE( z_esigwv3(jpi,jpj,jpk) )
557
558      z_gsigw3  = 0._wp   ;   z_gsigt3  = 0._wp   ;   z_gsi3w3  = 0._wp
559      z_esigt3  = 0._wp   ;   z_esigw3  = 0._wp 
560      z_esigtu3 = 0._wp   ;   z_esigtv3 = 0._wp   ;   z_esigtf3 = 0._wp
561      z_esigwu3 = 0._wp   ;   z_esigwv3 = 0._wp
562
563      DO ji = 1, jpi
564         DO jj = 1, jpj
565
566          IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
567             
568              zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
569                                                     ! could be changed by users but care must be taken to do so carefully
570              zzb = 1.0_wp-(zzb/hbatt(ji,jj))
571           
572              zzs = rn_zs / hbatt(ji,jj) 
573             
574              IF (rn_efold /= 0.0_wp) THEN
575                zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
576              ELSE
577                zsmth = 1.0_wp 
578              ENDIF
579               
580              DO jk = 1, jpk
581                z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
582                z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
583              ENDDO
584              z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth  )
585              z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth  )
586 
587          ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
588
589            DO jk = 1, jpk
590              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
591              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
592            END DO
593
594          ELSE  ! shallow water, z coordinates
595
596            DO jk = 1, jpk
597              z_gsigw3(ji,jj,jk) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
598              z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
599            END DO
600
601          ENDIF
602
603          DO jk = 1, jpk-1
604             z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
605             z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
606          END DO
607          z_esigw3(ji,jj,1  ) = 2.0_wp * (z_gsigt3(ji,jj,1  ) - z_gsigw3(ji,jj,1  ))
608          z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
609
610          ! Coefficients for vertical depth as the sum of e3w scale factors
611          z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
612          DO jk = 2, jpk
613             z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
614          END DO
615
616          DO jk = 1, jpk
617             gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
618             gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
619             gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
620          END DO
621
622        ENDDO   ! for all jj's
623      ENDDO    ! for all ji's
624
625      DO ji=1,jpi-1
626        DO jj=1,jpj-1
627
628          DO jk = 1, jpk
629                z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) / &
630                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
631                z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) / &
632                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
633                z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) +  &
634                                      hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / &
635                                    ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
636                z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) / &
637                                    ( hbatt(ji,jj)+hbatt(ji+1,jj) )
638                z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) / &
639                                    ( hbatt(ji,jj)+hbatt(ji,jj+1) )
640
641             e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
642             e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
643             e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
644             e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
645             !
646             e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
647             e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
648             e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
649          END DO
650
651        ENDDO
652      ENDDO
653      !
654      !                                               ! =============
655
656      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3)
657      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
658
659   END SUBROUTINE s_sf12
660
661   SUBROUTINE s_tanh()
662
663      !!----------------------------------------------------------------------
664      !!                  ***  ROUTINE s_tanh***
665      !!                     
666      !! ** Purpose :   stretch the s-coordinate system
667      !!
668      !! ** Method  :   s-coordinate stretch
669      !!
670      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
671      !!----------------------------------------------------------------------
672
673      INTEGER  ::   ji, jj, jk           ! dummy loop argument
674      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
675
676      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
677      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
678
679      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) )
680      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
681
682      z_gsigw  = 0._wp   ;   z_gsigt  = 0._wp   ;   z_gsi3w  = 0._wp
683      z_esigt  = 0._wp   ;   z_esigw  = 0._wp 
684
685      DO jk = 1, jpk
686        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
687        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
688      END DO
689      WRITE(*,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
690      !
691      ! Coefficients for vertical scale factors at w-, t- levels
692!!gm bug :  define it from analytical function, not like juste bellow....
693!!gm        or betteroffer the 2 possibilities....
694      DO jk = 1, jpk-1
695         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
696         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
697      END DO
698      z_esigw( 1 ) = 2._wp * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
699      z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
700      !
701      ! Coefficients for vertical depth as the sum of e3w scale factors
702      z_gsi3w(1) = 0.5_wp * z_esigw(1)
703      DO jk = 2, jpk
704         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
705      END DO
706!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
707      DO jk = 1, jpk
708         zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpk-1,wp)
709         zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpk-1,wp)
710         gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
711         gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
712         gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )
713      END DO
714!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
715      DO jj = 1, jpj
716         DO ji = 1, jpi
717            DO jk = 1, jpk
718              e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
719              e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
720              e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
721              e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) )
722              !
723              e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
724              e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
725              e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
726            END DO
727         END DO
728      END DO
729
730      DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w )
731      DEALLOCATE( z_esigt, z_esigw )
732
733   END SUBROUTINE s_tanh
734
735   FUNCTION fssig( pk ) RESULT( pf )
736      !!----------------------------------------------------------------------
737      !!                 ***  ROUTINE fssig ***
738      !!       
739      !! ** Purpose :   provide the analytical function in s-coordinate
740      !!         
741      !! ** Method  :   the function provide the non-dimensional position of
742      !!                T and W (i.e. between 0 and 1)
743      !!                T-points at integer values (between 1 and jpk)
744      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
745      !!----------------------------------------------------------------------
746      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
747      REAL(wp)             ::   pf   ! sigma value
748      !!----------------------------------------------------------------------
749      !
750      pf =   (   TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpk-1) + rn_thetb )  )   &
751         &     - TANH( rn_thetb * rn_theta                                )  )   &
752         & * (   COSH( rn_theta                           )                      &
753         &     + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) )  )              &
754         & / ( 2._wp * SINH( rn_theta ) )
755      !
756   END FUNCTION fssig
757
758
759   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
760      !!----------------------------------------------------------------------
761      !!                 ***  ROUTINE fssig1 ***
762      !!
763      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
764      !!
765      !! ** Method  :   the function provides the non-dimensional position of
766      !!                T and W (i.e. between 0 and 1)
767      !!                T-points at integer values (between 1 and jpk)
768      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
769      !!----------------------------------------------------------------------
770      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
771      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
772      REAL(wp)             ::   pf1   ! sigma value
773      !!----------------------------------------------------------------------
774      !
775      IF ( rn_theta == 0 ) then      ! uniform sigma
776         pf1 = - ( pk1 - 0.5_wp ) / REAL( jpk-1 )
777      ELSE                        ! stretched sigma
778         pf1 =   ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
779            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpk-1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta )  )  &
780            &        / ( 2._wp * TANH( 0.5_wp * rn_theta ) )  )
781      ENDIF
782      !
783   END FUNCTION fssig1
784
785
786   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
787      !!----------------------------------------------------------------------
788      !!                 ***  ROUTINE fgamma  ***
789      !!
790      !! ** Purpose :   provide analytical function for the s-coordinate
791      !!
792      !! ** Method  :   the function provides the non-dimensional position of
793      !!                T and W (i.e. between 0 and 1)
794      !!                T-points at integer values (between 1 and jpk)
795      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
796      !!
797      !!                This method allows the maintenance of fixed surface and or
798      !!                bottom cell resolutions (cf. geopotential coordinates)
799      !!                within an analytically derived stretched S-coordinate framework.
800      !!
801      !! Reference  :   Siddorn and Furner, in prep
802      !!----------------------------------------------------------------------
803      REAL(wp), INTENT(in   ) ::   pk1(jpk)       ! continuous "k" coordinate
804      REAL(wp)                ::   p_gamma(jpk)   ! stretched coordinate
805      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
806      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
807      REAL(wp), INTENT(in   ) ::   psmth       ! Smoothing parameter
808      REAL(wp)                ::   za1,za2,za3    ! local variables
809      REAL(wp)                ::   zn1,zn2        ! local variables
810      REAL(wp)                ::   za,zb,zx       ! local variables
811      integer                 ::   jk
812      !!----------------------------------------------------------------------
813      !
814
815      zn1  =  1./(jpk-1.)
816      zn2  =  1. -  zn1
817
818      za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 
819      za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
820      za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
821     
822      za = pzb - za3*(pzs-za1)-za2
823      za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
824      zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
825      zx = 1.0_wp-za/2.0_wp-zb
826 
827      DO jk = 1, jpk
828        p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp +  &
829                    & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
830                    &      (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
831        p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
832      ENDDO 
833
834      !
835   END FUNCTION fgamma
836
837END PROGRAM SCOORD_GEN
Note: See TracBrowser for help on using the repository browser.