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

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

1) Added tapering near equator if domain crosses equator
2) Added all variables that are needed for mesh_mask file to the output file.
File size: 35.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      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(*,*) 'scoord_gen : 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(*,*) '        Tapering in vicinity of equator              ln_eq_taper   = ',ln_eq_taper
69      WRITE(*,*) '        Horizontal Coordinate File                   cn_coord_hgr  = ',cn_coord_hgr
70      WRITE(*,*) '     Song and Haidvogel 1994 stretching              ln_s_sh94     = ',ln_s_sh94
71      WRITE(*,*) '        Song and Haidvogel 1994 stretching coefficients'
72      WRITE(*,*) '        surface control parameter (0<=rn_theta<=20)  rn_theta      = ',rn_theta
73      WRITE(*,*) '        bottom  control parameter (0<=rn_thetb<= 1)  rn_thetb      = ',rn_thetb
74      WRITE(*,*) '        stretching parameter (song and haidvogel)    rn_bb         = ',rn_bb
75      WRITE(*,*) '     Siddorn and Furner 2012 stretching              ln_s_sf12     = ',ln_s_sf12
76      WRITE(*,*) '        switching to sigma (T) or Z (F) at H<Hc      ln_sigcrit    = ',ln_sigcrit
77      WRITE(*,*) '        Siddorn and Furner 2012 stretching coefficients'
78      WRITE(*,*) '        stretchin parameter ( >1 surface; <1 bottom) rn_alpha      = ',rn_alpha
79      WRITE(*,*) '        e-fold length scale for transition region    rn_efold      = ',rn_efold
80      WRITE(*,*) '        Surface cell depth (Zs) (m)                  rn_zs         = ',rn_zs
81      WRITE(*,*) '        Bathymetry multiplier for Zb                 rn_zb_a       = ',rn_zb_a
82      WRITE(*,*) '        Offset for Zb                                rn_zb_b       = ',rn_zb_b
83      WRITE(*,*) '        Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
84
85     
86      ! Read in bathy, jpi and jpj from bathy.nc
87      CALL read_bathy()
88
89      !Allocate all other allocatable variables
90      ios = dom_oce_alloc()
91      IF( ios .NE. 0) THEN
92         WRITE(0,*) 'Unable to allocate all arrays'
93         STOP 1
94      ENDIF
95
96      hift(:,:) = rn_sbot_min                     ! set the minimum depth for the s-coordinate
97      hifu(:,:) = rn_sbot_min
98      hifv(:,:) = rn_sbot_min
99      hiff(:,:) = rn_sbot_min
100
101      !                                        ! set maximum ocean depth
102      bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
103
104      DO jj = 1, jpj
105         DO ji = 1, jpi
106           IF( bathy(ji,jj) > 0. )   bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
107         END DO
108      END DO
109      !                                        ! =============================
110      !                                        ! Define the envelop bathymetry   (hbatt)
111      !                                        ! =============================
112      ! use r-value to create hybrid coordinates
113      zenv(:,:) = bathy(:,:)
114      !
115      ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
116      DO jj = 1, jpj
117         DO ji = 1, jpi
118           IF( bathy(ji,jj) == 0. ) THEN
119             iip1 = MIN( ji+1, jpi )
120             ijp1 = MIN( jj+1, jpj )
121             iim1 = MAX( ji-1, 1 )
122             ijm1 = MAX( jj-1, 1 )
123             IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) +              &
124        &         bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0. ) THEN
125               zenv(ji,jj) = rn_sbot_min
126             ENDIF
127           ENDIF
128         END DO
129      END DO
130      !
131      ! smooth the bathymetry (if required)
132      scosrf(:,:) = 0.             ! ocean surface depth (here zero: no under ice-shelf sea)
133      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth
134      !
135      jl = 0
136      zrmax = 1.
137      !   
138      !     
139      ! set scaling factor used in reducing vertical gradients
140      zrfact = ( 1. - rn_rmax ) / ( 1. + rn_rmax )
141      !
142      ! initialise temporary evelope depth arrays
143      ztmpi1(:,:) = zenv(:,:)
144      ztmpi2(:,:) = zenv(:,:)
145      ztmpj1(:,:) = zenv(:,:)
146      ztmpj2(:,:) = zenv(:,:)
147      !
148      ! initialise temporary r-value arrays
149      zri(:,:) = 1.
150      zrj(:,:) = 1.
151
152      !                                                            ! ================ !
153      DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8 ) !  Iterative loop  !
154         !                                                         ! ================ !
155         jl = jl + 1
156         zrmax = 0.
157         ! we set zrmax from previous r-values (zri and zrj) first
158         ! if set after current r-value calculation (as previously)
159         ! we could exit DO WHILE prematurely before checking r-value
160         ! of current zenv
161         DO jj = 1, jpj 
162            DO ji = 1, jpi
163               zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
164            END DO
165         END DO
166         zri(:,:) = 0.
167         zrj(:,:) = 0.
168         DO jj = 1, jpj
169            DO ji = 1, jpi
170               iip1 = MIN(ji+1,jpi) 
171               ijp1 = MIN(jj+1,jpj)     
172               IF( (zenv(ji,jj) > 0.) .AND. (zenv(iip1,jj) > 0.)) THEN
173                  zri(ji,jj) = ( zenv(iip1,jj  ) - zenv(ji,jj) ) / ( zenv(iip1,jj  ) + zenv(ji,jj) )
174               END IF
175               IF( (zenv(ji,jj) > 0.) .AND. (zenv(ji,ijp1) > 0.)) THEN
176                  zrj(ji,jj) = ( zenv(ji  ,ijp1) - zenv(ji,jj) ) / ( zenv(ji  ,ijp1) + zenv(ji,jj) )
177               END IF
178               IF( zri(ji,jj) >  rn_rmax )   ztmpi1(ji  ,jj  ) = zenv(iip1,jj  ) * zrfact
179               IF( zri(ji,jj) < -rn_rmax )   ztmpi2(iip1,jj  ) = zenv(ji  ,jj  ) * zrfact
180               IF( zrj(ji,jj) >  rn_rmax )   ztmpj1(ji  ,jj  ) = zenv(ji  ,ijp1) * zrfact
181               IF( zrj(ji,jj) < -rn_rmax )   ztmpj2(ji  ,ijp1) = zenv(ji  ,jj  ) * zrfact
182            END DO
183         END DO
184         !
185         WRITE(*,*) 'zgr_sco :   iter= ',jl, ' rmax= ', zrmax
186         !
187         DO jj = 1, jpj
188            DO ji = 1, jpi
189               zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
190            END DO
191         END DO
192         !                                                  ! ================ !
193      END DO                                                !     End loop     !
194      !                                                     ! ================ !
195      DO jj = 1, jpj
196         DO ji = 1, jpi
197            zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
198         END DO
199      END DO
200      !
201      ! Envelope bathymetry saved in hbatt
202      hbatt(:,:) = zenv(:,:) 
203      ! TODO - get this section to work
204      IF( ln_eq_taper) THEN
205        CALL READ_GPHIT()
206        IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. ) THEN
207          CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
208          DO jj = 1, jpj
209             DO ji = 1, jpi
210                ztaper = EXP( -(gphit(ji,jj)/8.)**2. )
211                hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1. - ztaper )
212             END DO
213          END DO
214      ENDIF
215      !
216      !                                        ! ==============================
217      !                                        !   hbatu, hbatv, hbatf fields
218      !                                        ! ==============================
219      WRITE(*,*)
220      WRITE(*,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
221      hbatu(:,:) = rn_sbot_min
222      hbatv(:,:) = rn_sbot_min
223      hbatf(:,:) = rn_sbot_min
224      DO jj = 1, jpj-1
225        DO ji = 1, jpi-1   ! NO vector opt.
226           hbatu(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji+1,jj  ) )
227           hbatv(ji,jj) = 0.50 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1) )
228           hbatf(ji,jj) = 0.25 * ( hbatt(ji  ,jj) + hbatt(ji  ,jj+1)   &
229              &                     + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
230        END DO
231      END DO
232      !
233      zhbat(:,:) = hbatu(:,:)   
234      DO jj = 1, jpj
235         DO ji = 1, jpi
236            IF( hbatu(ji,jj) == 0. ) THEN
237               IF( zhbat(ji,jj) == 0. )   hbatu(ji,jj) = rn_sbot_min
238               IF( zhbat(ji,jj) /= 0. )   hbatu(ji,jj) = zhbat(ji,jj)
239            ENDIF
240         END DO
241      END DO
242      zhbat(:,:) = hbatv(:,:) 
243      DO jj = 1, jpj
244         DO ji = 1, jpi
245            IF( hbatv(ji,jj) == 0. ) THEN
246               IF( zhbat(ji,jj) == 0. )   hbatv(ji,jj) = rn_sbot_min
247               IF( zhbat(ji,jj) /= 0. )   hbatv(ji,jj) = zhbat(ji,jj)
248            ENDIF
249         END DO
250      END DO
251      zhbat(:,:) = hbatf(:,:)   
252      DO jj = 1, jpj
253         DO ji = 1, jpi
254            IF( hbatf(ji,jj) == 0. ) THEN
255               IF( zhbat(ji,jj) == 0. )   hbatf(ji,jj) = rn_sbot_min
256               IF( zhbat(ji,jj) /= 0. )   hbatf(ji,jj) = zhbat(ji,jj)
257            ENDIF
258         END DO
259      END DO
260
261!!bug:  key_helsinki a verifer
262      hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
263      hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
264      hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
265      hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
266
267         WRITE(*,*) ' MAX val hif   t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ),  &
268            &                        ' u ',   MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
269         WRITE(*,*) ' MIN val hif   t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ),  &
270            &                        ' u ',   MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
271         WRITE(*,*) ' MAX val hbat  t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ),  &
272            &                        ' u ',   MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
273         WRITE(*,*) ' MIN val hbat  t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ),  &
274            &                        ' u ',   MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
275     
276      ! Create the output file
277      CALL make_coord_file()
278
279      !                                            ! =======================
280      !                                            !   s-coordinate fields     (gdep., e3.)
281      !                                            ! =======================
282      !
283      ! non-dimensional "sigma" for model level depth at w- and t-levels
284
285
286!========================================================================
287! Song and Haidvogel  1994 (ln_s_sh94=T)
288! Siddorn and Furner 2012 (ln_sf12=T)
289! or  tanh function       (both false)                   
290! To reduce memory loop over jpk and write each level to file
291!========================================================================
292      IF      ( ln_s_sh94 ) THEN
293                           CALL s_sh94()
294      ELSE IF ( ln_s_sf12 ) THEN
295                           CALL s_sf12()
296      ELSE                 
297                           CALL s_tanh()
298      ENDIF
299
300      CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) )
301      CALL check_nf90( nf90_close(ncout) )
302
303 
304      !
305END PROGRAM SCOORD_GEN
306
307!!======================================================================
308   SUBROUTINE s_sh94()
309
310      !!----------------------------------------------------------------------
311      !!                  ***  ROUTINE s_sh94  ***
312      !!                     
313      !! ** Purpose :   stretch the s-coordinate system
314      !!
315      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
316      !!                mixed S/sigma coordinate
317      !!
318      !! Reference : Song and Haidvogel 1994.
319      !!----------------------------------------------------------------------
320      !
321      USE utils
322      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
323      !
324      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
325      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
326      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
327      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
328
329      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
330      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
331      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
332      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
333      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
334 
335      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
336      z_esigt3  = 0.   ;   z_esigw3  = 0. 
337      z_esigt3p1  = 0. ;   z_esigw3p1  = 0. 
338      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
339      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
340
341      DO jk = 1,jpk
342         DO ji = 1, jpi
343            DO jj = 1, jpj
344
345               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
346                     z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb )
347                     z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb )
348                     z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp)    , rn_bb )
349               ELSE ! shallow water, uniform sigma
350                     z_gsigw3(ji,jj) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
351                     z_gsigw3p1(ji,jj) =   REAL(jk,wp)            / REAL(jpk-1,wp)
352                     z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp)
353               ENDIF
354               !
355             !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
356               IF( jk .EQ. 1) THEN
357                  z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) )
358                  z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
359               ELSE
360                  z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
361                  z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
362               ENDIF
363               IF( jk .EQ. jpk) THEN
364                  z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) )
365               ELSE
366                  z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
367               ENDIF
368               !
369
370               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
371               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
372               gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft )
373               gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw )
374               gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft )
375              !
376            END DO   ! for all jj's
377         END DO    ! for all ji's
378
379         DO ji = 1, jpi-1
380            DO jj = 1, jpj-1
381               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) )   &
382                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
383               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) )   &
384                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
385               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj)     &
386                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) )   &
387                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
388               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) )   &
389                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
390               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) )   &
391                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
392               !
393               e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
394               e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
395               e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
396               e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
397               !
398               e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
399               e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
400               e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
401            END DO
402         END DO
403
404         z_gsigt3m1 = z_gsigt3
405         z_gsiw3m1 = z_gsiw3
406 
407         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
408         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
409         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
410         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
411         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
412         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
413         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
414       
415         CALL write_netcdf_vars(jk)
416         DO jj = 1, jpj
417            DO ji = 1, jpi
418                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
419                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
420            END DO
421         END DO
422      END DO !End of loop over jk
423
424      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
425      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
426
427   END SUBROUTINE s_sh94
428
429   SUBROUTINE s_sf12
430
431      !!----------------------------------------------------------------------
432      !!                  ***  ROUTINE s_sf12 ***
433      !!                     
434      !! ** Purpose :   stretch the s-coordinate system
435      !!
436      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
437      !!                mixed S/sigma/Z coordinate
438      !!
439      !!                This method allows the maintenance of fixed surface and or
440      !!                bottom cell resolutions (cf. geopotential coordinates)
441      !!                within an analytically derived stretched S-coordinate framework.
442      !!
443      !!
444      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
445      !!----------------------------------------------------------------------
446      !
447      USE utils
448      REAL(wp) ::   zsmth               ! smoothing around critical depth
449      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
450      !
451      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
452      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
453      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
454      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
455
456      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
457      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
458      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
459      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
460      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
461
462      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
463      z_esigt3  = 0.   ;   z_esigw3  = 0. 
464      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
465      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
466     
467
468      DO jk = 1,jpk
469         DO ji = 1, jpi
470            DO jj = 1, jpj
471             IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
472             
473                zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
474                                                     ! could be changed by users but care must be taken to do so carefully
475                zzb = 1.0-(zzb/hbatt(ji,jj))
476           
477                zzs = rn_zs / hbatt(ji,jj) 
478             
479                IF (rn_efold /= 0.0) THEN
480                   zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
481                ELSE
482                   zsmth = 1.0 
483                ENDIF
484                 
485                z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
486                z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)
487                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
488                z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth  )
489                z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(ji,jj), zzb, zzs, zsmth  )
490                z_gsigt3(ji,jj) = fgamma( z_gsigt3(ji,jj), zzb, zzs, zsmth  )
491 
492             ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma
493
494                z_gsigw3(ji,jj) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
495                z_gsigw3p1(ji,jj) =  REAL(jk,wp)     /REAL(jpk-1,wp)
496                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
497
498             ELSE  ! shallow water, z coordinates
499
500               z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
501               z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
502               z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
503
504             ENDIF
505
506             !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
507             IF( jk .EQ. 1) THEN
508                z_esigw3(ji,jj  ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ))
509                z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
510             ELSE
511                z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
512                z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
513             ENDIF
514             IF( jk .EQ. jpk) THEN
515                z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj))
516             ELSE
517                z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
518             ENDIF
519
520             gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj)
521             gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj)
522             gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj)
523
524            ENDDO   ! for all jj's
525         ENDDO    ! for all ji's
526
527         DO ji=1,jpi-1
528            DO jj=1,jpj-1
529
530               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / &
531                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
532               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / &
533                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
534               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) +  &
535                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / &
536                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
537               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / &
538                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
539               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / &
540                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
541
542               e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj)
543               e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj)
544               e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj)
545               e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj)
546               !
547               e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj)
548               e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj)
549               e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj)
550            ENDDO
551         ENDDO
552         ! Keep some arrays for next level
553         z_gsigt3m1 = z_gsigt3
554         z_gsiw3m1 = z_gsiw3
555 
556         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
557         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
558         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
559         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
560         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
561         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
562         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
563         
564         CALL write_netcdf_vars(jk)
565
566         DO jj = 1, jpj
567            DO ji = 1, jpi
568                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
569                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
570            END DO
571         END DO
572
573      ENDDO ! End of loop over jk
574
575      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
576      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
577
578   END SUBROUTINE s_sf12
579
580   SUBROUTINE s_tanh()
581
582      !!----------------------------------------------------------------------
583      !!                  ***  ROUTINE s_tanh***
584      !!                     
585      !! ** Purpose :   stretch the s-coordinate system
586      !!
587      !! ** Method  :   s-coordinate stretch
588      !!
589      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
590      !!----------------------------------------------------------------------
591
592      USE utils
593      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
594
595      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
596      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
597
598      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) )
599      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
600
601      z_gsigw  = 0.   ;   z_gsigt  = 0.   ;   z_gsi3w  = 0.
602      z_esigt  = 0.   ;   z_esigw  = 0. 
603
604      DO jk = 1, jpk
605        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 )
606        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
607      END DO
608      WRITE(*,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
609      !
610      ! Coefficients for vertical scale factors at w-, t- levels
611!!gm bug :  define it from analytical function, not like juste bellow....
612!!gm        or betteroffer the 2 possibilities....
613      DO jk = 1, jpk-1
614         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
615         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
616      END DO
617      z_esigw( 1 ) = 2. * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
618      z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) )
619      !
620      ! Coefficients for vertical depth as the sum of e3w scale factors
621      z_gsi3w(1) = 0.5 * z_esigw(1)
622      DO jk = 2, jpk
623         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
624      END DO
625!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
626      DO jk = 1, jpk
627      END DO
628!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
629      DO jk = 1, jpk
630         DO jj = 1, jpj
631            DO ji = 1, jpi
632              zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
633              zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
634              gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft )
635              gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw )
636              gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft )
637              e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
638              e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
639              e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
640              e3f_0(ji,jj) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) )
641              !
642              e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
643              e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
644              e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
645              IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
646              IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
647            END DO
648         END DO
649         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
650         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
651         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
652         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
653         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
654         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
655         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
656 
657         CALL write_netcdf_vars(jk)
658      ENDDO ! End of loop over jk
659
660
661      DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w )
662      DEALLOCATE( z_esigt, z_esigw )
663
664   END SUBROUTINE s_tanh
665
666   FUNCTION fssig( pk ) RESULT( pf )
667      !!----------------------------------------------------------------------
668      !!                 ***  ROUTINE fssig ***
669      !!       
670      !! ** Purpose :   provide the analytical function in s-coordinate
671      !!         
672      !! ** Method  :   the function provide the non-dimensional position of
673      !!                T and W (i.e. between 0 and 1)
674      !!                T-points at integer values (between 1 and jpk)
675      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
676      !!----------------------------------------------------------------------
677      USE utils, ONLY : wp,rn_theta,rn_thetb,jpk
678      IMPLICIT NONE
679      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
680      REAL(wp)             ::   pf   ! sigma value
681      !!----------------------------------------------------------------------
682      !
683      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb )  )   &
684         &     - TANH( rn_thetb * rn_theta                                )  )   &
685         & * (   COSH( rn_theta                           )                      &
686         &     + COSH( rn_theta * ( 2. * rn_thetb - 1. ) )  )              &
687         & / ( 2. * SINH( rn_theta ) )
688      !
689   END FUNCTION fssig
690
691
692   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
693      !!----------------------------------------------------------------------
694      !!                 ***  ROUTINE fssig1 ***
695      !!
696      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
697      !!
698      !! ** Method  :   the function provides the non-dimensional position of
699      !!                T and W (i.e. between 0 and 1)
700      !!                T-points at integer values (between 1 and jpk)
701      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
702      !!----------------------------------------------------------------------
703      USE utils, ONLY : wp, jpk, rn_theta
704      IMPLICIT NONE
705      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
706      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
707      REAL(wp)             ::   pf1   ! sigma value
708      !!----------------------------------------------------------------------
709      !
710      IF ( rn_theta == 0 ) then      ! uniform sigma
711         pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 )
712      ELSE                        ! stretched sigma
713         pf1 =   ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
714            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta )  )  &
715            &        / ( 2. * TANH( 0.5 * rn_theta ) )  )
716      ENDIF
717      !
718   END FUNCTION fssig1
719
720
721   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
722      !!----------------------------------------------------------------------
723      !!                 ***  ROUTINE fgamma  ***
724      !!
725      !! ** Purpose :   provide analytical function for the s-coordinate
726      !!
727      !! ** Method  :   the function provides the non-dimensional position of
728      !!                T and W (i.e. between 0 and 1)
729      !!                T-points at integer values (between 1 and jpk)
730      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
731      !!
732      !!                This method allows the maintenance of fixed surface and or
733      !!                bottom cell resolutions (cf. geopotential coordinates)
734      !!                within an analytically derived stretched S-coordinate framework.
735      !!
736      !! Reference  :   Siddorn and Furner, in prep
737      !!----------------------------------------------------------------------
738      USE utils, ONLY : jpk,wp,rn_alpha
739      IMPLICIT NONE
740      REAL(wp), INTENT(in   ) ::   pk1           ! continuous "k" coordinate
741      REAL(wp)                ::   p_gamma       ! stretched coordinate
742      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
743      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
744      REAL(wp), INTENT(in   ) ::   psmth         ! Smoothing parameter
745      REAL(wp)                ::   za1,za2,za3   ! local variables
746      REAL(wp)                ::   zn1,zn2       ! local variables
747      REAL(wp)                ::   za,zb,zx      ! local variables
748      !!----------------------------------------------------------------------
749      !
750
751      zn1  =  1./(jpk-1.)
752      zn2  =  1. -  zn1
753
754      za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) 
755      za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0)
756      za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1)
757     
758      za = pzb - za3*(pzs-za1)-za2
759      za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) )
760      zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1)
761      zx = 1.0-za/2.0-zb
762
763      p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 +  &
764                  & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- &
765                  &      (rn_alpha+1.0)*pk1**(rn_alpha+2.0) )
766      p_gamma = p_gamma*psmth+pk1*(1.0-psmth)
767
768      !
769   END FUNCTION fgamma
770
Note: See TracBrowser for help on using the repository browser.