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

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

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