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

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

Removed use of 3D variables.
Compiles and runs with no errors but still need to check output

File size: 36.5 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     
289      ! Create the output file
290      CALL make_coord_file()
291
292      !                                            ! =======================
293      !                                            !   s-coordinate fields     (gdep., e3.)
294      !                                            ! =======================
295      !
296      ! non-dimensional "sigma" for model level depth at w- and t-levels
297
298
299!========================================================================
300! Song and Haidvogel  1994 (ln_s_sh94=T)
301! Siddorn and Furner 2012 (ln_sf12=T)
302! or  tanh function       (both false)                   
303! To reduce memory loop over jpk and write each level to file
304!========================================================================
305      IF      ( ln_s_sh94 ) THEN
306                           CALL s_sh94()
307      ELSE IF ( ln_s_sf12 ) THEN
308                           CALL s_sf12()
309      ELSE                 
310                           CALL s_tanh()
311      ENDIF
312
313      CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy) )
314      CALL check_nf90( nf90_close(ncout) )
315
316 
317      !
318END PROGRAM SCOORD_GEN
319
320!!======================================================================
321   SUBROUTINE s_sh94()
322
323      !!----------------------------------------------------------------------
324      !!                  ***  ROUTINE s_sh94  ***
325      !!                     
326      !! ** Purpose :   stretch the s-coordinate system
327      !!
328      !! ** Method  :   s-coordinate stretch using the Song and Haidvogel 1994
329      !!                mixed S/sigma coordinate
330      !!
331      !! Reference : Song and Haidvogel 1994.
332      !!----------------------------------------------------------------------
333      !
334      USE utils
335      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
336      !
337      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
338      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
339      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
340      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
341
342      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
343      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
344      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
345      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
346      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
347 
348      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
349      z_esigt3  = 0.   ;   z_esigw3  = 0. 
350      z_esigt3p1  = 0. ;   z_esigw3p1  = 0. 
351      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
352      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
353
354      DO jk = 1,jpk
355         DO ji = 1, jpi
356            DO jj = 1, jpj
357
358               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
359                     z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb )
360                     z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb )
361                     z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp)    , rn_bb )
362               ELSE ! shallow water, uniform sigma
363                     z_gsigw3(ji,jj) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
364                     z_gsigw3p1(ji,jj) =   REAL(jk,wp)            / REAL(jpk-1,wp)
365                     z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp)
366               ENDIF
367               !
368             !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
369               IF( jk .EQ. 1) THEN
370                  z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) )
371                  z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
372               ELSE
373                  z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
374                  z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
375               ENDIF
376               IF( jk .EQ. jpk) THEN
377                  z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) )
378               ELSE
379                  z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
380               ENDIF
381               !
382
383               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
384               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
385               gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft )
386               gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw )
387               gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft )
388              !
389            END DO   ! for all jj's
390         END DO    ! for all ji's
391
392         DO ji = 1, jpim1
393            DO jj = 1, jpjm1
394               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) )   &
395                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
396               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) )   &
397                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
398               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj)     &
399                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) )   &
400                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
401               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) )   &
402                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
403               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) )   &
404                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
405               !
406               e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
407               e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
408               e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
409               e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
410               !
411               e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
412               e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
413               e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
414            END DO
415         END DO
416
417         z_gsigt3m1 = z_gsigt3
418         z_gsiw3m1 = z_gsiw3
419 
420         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
421         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
422         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
423         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
424         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
425         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
426         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
427       
428         CALL write_netcdf_vars(jk)
429         DO jj = 1, jpj
430            DO ji = 1, jpi
431                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
432                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
433            END DO
434         END DO
435      END DO !End of loop over jk
436
437      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
438      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
439
440   END SUBROUTINE s_sh94
441
442   SUBROUTINE s_sf12
443
444      !!----------------------------------------------------------------------
445      !!                  ***  ROUTINE s_sf12 ***
446      !!                     
447      !! ** Purpose :   stretch the s-coordinate system
448      !!
449      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
450      !!                mixed S/sigma/Z coordinate
451      !!
452      !!                This method allows the maintenance of fixed surface and or
453      !!                bottom cell resolutions (cf. geopotential coordinates)
454      !!                within an analytically derived stretched S-coordinate framework.
455      !!
456      !!
457      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
458      !!----------------------------------------------------------------------
459      !
460      USE utils
461      REAL(wp) ::   zsmth               ! smoothing around critical depth
462      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
463      !
464      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
465      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
466      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
467      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
468
469      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
470      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
471      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
472      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
473      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
474
475      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
476      z_esigt3  = 0.   ;   z_esigw3  = 0. 
477      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
478      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
479     
480
481      DO jk = 1,jpk
482         DO ji = 1, jpi
483            DO jj = 1, jpj
484             IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
485             
486                zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
487                                                     ! could be changed by users but care must be taken to do so carefully
488                zzb = 1.0-(zzb/hbatt(ji,jj))
489           
490                zzs = rn_zs / hbatt(ji,jj) 
491             
492                IF (rn_efold /= 0.0) THEN
493                   zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
494                ELSE
495                   zsmth = 1.0 
496                ENDIF
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                z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth  )
502                z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(ji,jj), zzb, zzs, zsmth  )
503                z_gsigt3(ji,jj) = fgamma( z_gsigt3(ji,jj), zzb, zzs, zsmth  )
504 
505             ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma
506
507                z_gsigw3(ji,jj) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
508                z_gsigw3p1(ji,jj) =  REAL(jk,wp)     /REAL(jpk-1,wp)
509                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
510
511             ELSE  ! shallow water, z coordinates
512
513               z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
514               z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
515               z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
516
517             ENDIF
518
519             !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
520             IF( jk .EQ. 1) THEN
521                z_esigw3(ji,jj  ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ))
522                z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
523             ELSE
524                z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
525                z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
526             ENDIF
527             IF( jk .EQ. jpk) THEN
528                z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj))
529             ELSE
530                z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
531             ENDIF
532
533             gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj)
534             gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj)
535             gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj)
536
537            ENDDO   ! for all jj's
538         ENDDO    ! for all ji's
539
540         DO ji=1,jpi-1
541            DO jj=1,jpj-1
542
543               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / &
544                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
545               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / &
546                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
547               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) +  &
548                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / &
549                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
550               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / &
551                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
552               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / &
553                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
554
555               e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj)
556               e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj)
557               e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj)
558               e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj)
559               !
560               e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj)
561               e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj)
562               e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj)
563            ENDDO
564         ENDDO
565         ! Keep some arrays for next level
566         z_gsigt3m1 = z_gsigt3
567         z_gsiw3m1 = z_gsiw3
568 
569         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
570         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
571         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
572         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
573         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
574         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
575         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
576         
577         WRITE(*,*) 'Writing level ',jk,' to file'
578         CALL write_netcdf_vars(jk)
579         WRITE(*,*) 'Written level ',jk,' to file'
580      ENDDO ! End of loop over jk
581
582      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
583      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
584
585   END SUBROUTINE s_sf12
586
587   SUBROUTINE s_tanh()
588
589      !!----------------------------------------------------------------------
590      !!                  ***  ROUTINE s_tanh***
591      !!                     
592      !! ** Purpose :   stretch the s-coordinate system
593      !!
594      !! ** Method  :   s-coordinate stretch
595      !!
596      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
597      !!----------------------------------------------------------------------
598
599      USE utils
600      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
601
602      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
603      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
604
605      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) )
606      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
607
608      z_gsigw  = 0.   ;   z_gsigt  = 0.   ;   z_gsi3w  = 0.
609      z_esigt  = 0.   ;   z_esigw  = 0. 
610
611      DO jk = 1, jpk
612        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 )
613        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
614      END DO
615      WRITE(*,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
616      !
617      ! Coefficients for vertical scale factors at w-, t- levels
618!!gm bug :  define it from analytical function, not like juste bellow....
619!!gm        or betteroffer the 2 possibilities....
620      DO jk = 1, jpk-1
621         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
622         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
623      END DO
624      z_esigw( 1 ) = 2. * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
625      z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) )
626      !
627      ! Coefficients for vertical depth as the sum of e3w scale factors
628      z_gsi3w(1) = 0.5 * z_esigw(1)
629      DO jk = 2, jpk
630         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
631      END DO
632!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
633      DO jk = 1, jpk
634      END DO
635!!gm: e3uw, e3vw can be suppressed  (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
636      DO jk = 1, jpk
637         DO jj = 1, jpj
638            DO ji = 1, jpi
639              zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
640              zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
641              gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigt(jk) + hift(ji,jj)*zcoeft )
642              gdepw_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsigw(jk) + hift(ji,jj)*zcoefw )
643              gdep3w_0(:,:) = ( scosrf(ji,jj) + (hbatt(ji,jj)-hift(ji,jj))*z_gsi3w(jk) + hift(ji,jj)*zcoeft )
644              e3t_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
645              e3u_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
646              e3v_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
647              e3f_0(ji,jj) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpk-1,wp) )
648              !
649              e3w_0(ji,jj) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpk-1,wp) )
650              e3uw_0(ji,jj) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpk-1,wp) )
651              e3vw_0(ji,jj) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpk-1,wp) )
652              IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
653              IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
654            END DO
655         END DO
656         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
657         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
658         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
659         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
660         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
661         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
662         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
663 
664         CALL write_netcdf_vars(jk)
665      ENDDO ! End of loop over jk
666
667
668      DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w )
669      DEALLOCATE( z_esigt, z_esigw )
670
671   END SUBROUTINE s_tanh
672
673   FUNCTION fssig( pk ) RESULT( pf )
674      !!----------------------------------------------------------------------
675      !!                 ***  ROUTINE fssig ***
676      !!       
677      !! ** Purpose :   provide the analytical function in s-coordinate
678      !!         
679      !! ** Method  :   the function provide the non-dimensional position of
680      !!                T and W (i.e. between 0 and 1)
681      !!                T-points at integer values (between 1 and jpk)
682      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
683      !!----------------------------------------------------------------------
684      USE utils, ONLY : wp
685      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
686      REAL(wp)             ::   pf   ! sigma value
687      !!----------------------------------------------------------------------
688      !
689      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb )  )   &
690         &     - TANH( rn_thetb * rn_theta                                )  )   &
691         & * (   COSH( rn_theta                           )                      &
692         &     + COSH( rn_theta * ( 2. * rn_thetb - 1. ) )  )              &
693         & / ( 2. * SINH( rn_theta ) )
694      !
695   END FUNCTION fssig
696
697
698   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
699      !!----------------------------------------------------------------------
700      !!                 ***  ROUTINE fssig1 ***
701      !!
702      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
703      !!
704      !! ** Method  :   the function provides the non-dimensional position of
705      !!                T and W (i.e. between 0 and 1)
706      !!                T-points at integer values (between 1 and jpk)
707      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
708      !!----------------------------------------------------------------------
709      USE utils, ONLY : wp
710      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
711      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
712      REAL(wp)             ::   pf1   ! sigma value
713      !!----------------------------------------------------------------------
714      !
715      IF ( rn_theta == 0 ) then      ! uniform sigma
716         pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 )
717      ELSE                        ! stretched sigma
718         pf1 =   ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
719            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta )  )  &
720            &        / ( 2. * TANH( 0.5 * rn_theta ) )  )
721      ENDIF
722      !
723   END FUNCTION fssig1
724
725
726   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
727      !!----------------------------------------------------------------------
728      !!                 ***  ROUTINE fgamma  ***
729      !!
730      !! ** Purpose :   provide analytical function for the s-coordinate
731      !!
732      !! ** Method  :   the function provides the non-dimensional position of
733      !!                T and W (i.e. between 0 and 1)
734      !!                T-points at integer values (between 1 and jpk)
735      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
736      !!
737      !!                This method allows the maintenance of fixed surface and or
738      !!                bottom cell resolutions (cf. geopotential coordinates)
739      !!                within an analytically derived stretched S-coordinate framework.
740      !!
741      !! Reference  :   Siddorn and Furner, in prep
742      !!----------------------------------------------------------------------
743      USE utils, ONLY : jpk,jk,wp
744      REAL(wp), INTENT(in   ) ::   pk1           ! continuous "k" coordinate
745      REAL(wp)                ::   p_gamma       ! stretched coordinate
746      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
747      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
748      REAL(wp), INTENT(in   ) ::   psmth         ! Smoothing parameter
749      REAL(wp)                ::   za1,za2,za3   ! local variables
750      REAL(wp)                ::   zn1,zn2       ! local variables
751      REAL(wp)                ::   za,zb,zx      ! local variables
752      !!----------------------------------------------------------------------
753      !
754
755      zn1  =  1./(jpk-1.)
756      zn2  =  1. -  zn1
757
758      za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) 
759      za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0)
760      za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1)
761     
762      za = pzb - za3*(pzs-za1)-za2
763      za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) )
764      zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1)
765      zx = 1.0-za/2.0-zb
766 
767      p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 +  &
768                  & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- &
769                  &      (rn_alpha+1.0)*pk1**(rn_alpha+2.0) )
770      p_gamma = p_gamma*psmth+pk1*(1.0-psmth)
771
772      !
773   END FUNCTION fgamma
774
Note: See TracBrowser for help on using the repository browser.