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

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

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