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 utils/tools_AGRIF_CMEMS_2020/SCOORD_GEN/src – NEMO

source: utils/tools_AGRIF_CMEMS_2020/SCOORD_GEN/src/scoord_gen.F90 @ 10978

Last change on this file since 10978 was 5993, checked in by timgraham, 8 years ago

Applied changes as suggested by reviewer
File size: 30.8 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_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
342      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
343
344      DO jk = 1,jpk
345         DO ji = 1, jpi
346            DO jj = 1, jpj
347
348               IF( hbatt(ji,jj) > rn_hc ) THEN    !deep water, stretched sigma
349                     z_gsigw3(ji,jj) = -fssig1( REAL(jk,wp)-0.5, rn_bb )
350                     z_gsigw3p1(ji,jj) = -fssig1( REAL(jk+1,wp)-0.5, rn_bb )
351                     z_gsigt3(ji,jj) = -fssig1( REAL(jk,wp)    , rn_bb )
352               ELSE ! shallow water, uniform sigma
353                     z_gsigw3(ji,jj) =   REAL(jk-1,wp)            / REAL(jpk-1,wp)
354                     z_gsigw3p1(ji,jj) =   REAL(jk,wp)            / REAL(jpk-1,wp)
355                     z_gsigt3(ji,jj) = ( REAL(jk-1,wp) + 0.5 ) / REAL(jpk-1,wp)
356               ENDIF
357               !
358             !gsi3w3m1 & gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
359               IF( jk .EQ. 1) THEN
360                  z_esigw3(ji,jj ) = 2. * ( z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ) )
361                  z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj)
362               ELSE
363                  z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
364                  z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
365               ENDIF
366               IF( jk .EQ. jpk) THEN
367                  z_esigt3(ji,jj) = 2. * ( z_gsigt3(ji,jj) - z_gsigw3(ji,jj) )
368               ELSE
369                  z_esigt3(ji,jj ) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
370               ENDIF
371               !
372
373               zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpk-1,wp)
374               zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpk-1,wp)
375               gdept_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj)+rn_hc*zcoeft )
376               gdepw_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj)+rn_hc*zcoefw )
377               gdep3w_0(ji,jj) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj)+rn_hc*zcoeft )
378              !
379            END DO   ! for all jj's
380         END DO    ! for all ji's
381
382         DO ji = 1, jpi-1
383            DO jj = 1, jpj-1
384               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) )   &
385                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
386               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) )   &
387                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
388               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj)     &
389                  &                + hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) )   &
390                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
391               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) )   &
392                  &              / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
393               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) )   &
394                  &              / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
395               !
396               e3t_0(ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
397               e3u_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
398               e3v_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
399               e3f_0(ji,jj) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
400               !
401               e3w_0 (ji,jj) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj) + rn_hc/REAL(jpk-1,wp) )
402               e3uw_0(ji,jj) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
403               e3vw_0(ji,jj) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj) + rn_hc/REAL(jpk-1,wp) )
404            END DO
405         END DO
406
407         z_gsigt3m1 = z_gsigt3
408         z_gsi3w3m1 = z_gsi3w3
409 
410         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
411         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
412         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
413         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
414         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
415         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
416         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
417       
418         CALL write_netcdf_3d_vars(jk)
419         DO jj = 1, jpj
420            DO ji = 1, jpi
421                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
422                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
423            END DO
424         END DO
425      END DO !End of loop over jk
426
427      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
428      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
429
430   END SUBROUTINE s_sh94
431
432   SUBROUTINE s_sf12
433
434      !!----------------------------------------------------------------------
435      !!                  ***  ROUTINE s_sf12 ***
436      !!                     
437      !! ** Purpose :   stretch the s-coordinate system
438      !!
439      !! ** Method  :   s-coordinate stretch using the Siddorn and Furner 2012?
440      !!                mixed S/sigma/Z coordinate
441      !!
442      !!                This method allows the maintenance of fixed surface and or
443      !!                bottom cell resolutions (cf. geopotential coordinates)
444      !!                within an analytically derived stretched S-coordinate framework.
445      !!
446      !!
447      !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
448      !!----------------------------------------------------------------------
449      !
450      USE utils
451      REAL(wp) ::   zsmth               ! smoothing around critical depth
452      REAL(wp) ::   zzs, zzb           ! Surface and bottom cell thickness in sigma space
453      !
454      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
455      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1
456      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigt3, z_esigw3, z_esigtu3 
457      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3           
458
459      ALLOCATE( z_gsigw3(jpi,jpj), z_gsigt3(jpi,jpj), z_gsi3w3(jpi,jpj) )
460      ALLOCATE( z_esigt3(jpi,jpj), z_esigw3(jpi,jpj), z_esigtu3(jpi,jpj)) 
461      ALLOCATE( z_esigtv3(jpi,jpj), z_esigtf3(jpi,jpj), z_esigwu3(jpi,jpj)) 
462      ALLOCATE( z_esigwv3(jpi,jpj), z_gsigw3p1(jpi,jpj), z_gsigt3m1(jpi,jpj) )
463      ALLOCATE( z_gsi3w3m1(jpi,jpj) )
464
465      z_gsigw3  = 0.   ;   z_gsigt3  = 0.   ;   z_gsi3w3  = 0.
466      z_esigt3  = 0.   ;   z_esigw3  = 0. 
467      z_esigtu3 = 0.   ;   z_esigtv3 = 0.   ;   z_esigtf3 = 0.
468      z_esigwu3 = 0.   ;   z_esigwv3 = 0.
469     
470
471      DO jk = 1,jpk
472         DO ji = 1, jpi
473            DO jj = 1, jpj
474             IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
475             
476                zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b   ! this forces a linear bottom cell depth relationship with H,.
477                                                     ! could be changed by users but care must be taken to do so carefully
478                zzb = 1.0-(zzb/hbatt(ji,jj))
479           
480                zzs = rn_zs / hbatt(ji,jj) 
481             
482                IF (rn_efold /= 0.0) THEN
483                   zsmth   = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
484                ELSE
485                   zsmth = 1.0 
486                ENDIF
487                 
488                z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)
489                z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)
490                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
491                z_gsigw3(ji,jj) = fgamma( z_gsigw3(ji,jj), zzb, zzs, zsmth  )
492                z_gsigw3p1(ji,jj) = fgamma( z_gsigw3p1(ji,jj), zzb, zzs, zsmth  )
493                z_gsigt3(ji,jj) = fgamma( z_gsigt3(ji,jj), zzb, zzs, zsmth  )
494 
495             ELSE IF(ln_sigcrit) THEN ! shallow water, uniform sigma
496
497                z_gsigw3(ji,jj) =  REAL(jk-1,wp)     /REAL(jpk-1,wp)
498                z_gsigw3p1(ji,jj) =  REAL(jk,wp)     /REAL(jpk-1,wp)
499                z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
500
501             ELSE  ! shallow water, z coordinates
502
503               z_gsigw3(ji,jj) =  REAL(jk-1,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
504               z_gsigw3p1(ji,jj) =  REAL(jk,wp)        /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
505               z_gsigt3(ji,jj) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
506
507             ENDIF
508
509             !gsi3w3m1 & z_gsigt3m1 only used if jk /= 1 and is set at the end of the loop over jk
510             IF( jk .EQ. 1) THEN
511                z_esigw3(ji,jj  ) = 2.0 * (z_gsigt3(ji,jj ) - z_gsigw3(ji,jj ))
512                z_gsi3w3(ji,jj) = 0.5 * z_esigw3(ji,jj) 
513             ELSE
514                z_esigw3(ji,jj) = z_gsigt3(ji,jj) - z_gsigt3m1(ji,jj)
515                z_gsi3w3(ji,jj) = z_gsi3w3m1(ji,jj) + z_esigw3(ji,jj)
516             ENDIF
517             IF( jk .EQ. jpk) THEN
518                z_esigt3(ji,jj) = 2.0 * (z_gsigt3(ji,jj) - z_gsigw3(ji,jj))
519             ELSE
520                z_esigt3(ji,jj) = z_gsigw3p1(ji,jj) - z_gsigw3(ji,jj)
521             ENDIF
522
523             gdept_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj)
524             gdepw_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj)
525             gdep3w_0(ji,jj) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj)
526
527            ENDDO   ! for all jj's
528         ENDDO    ! for all ji's
529
530         DO ji=1,jpi-1
531            DO jj=1,jpj-1
532
533               z_esigtu3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) ) / &
534                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
535               z_esigtv3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1) ) / &
536                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
537               z_esigtf3(ji,jj) = ( hbatt(ji,jj)*z_esigt3(ji,jj)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj) +  &
538                                       hbatt(ji,jj+1)*z_esigt3(ji,jj+1)+hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1) ) / &
539                                     ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
540               z_esigwu3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj) ) / &
541                                     ( hbatt(ji,jj)+hbatt(ji+1,jj) )
542               z_esigwv3(ji,jj) = ( hbatt(ji,jj)*z_esigw3(ji,jj)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1) ) / &
543                                     ( hbatt(ji,jj)+hbatt(ji,jj+1) )
544
545               e3t_0(ji,jj)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj)
546               e3u_0(ji,jj)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj)
547               e3v_0(ji,jj)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj)
548               e3f_0(ji,jj)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj)
549               !
550               e3w_0(ji,jj)=hbatt(ji,jj)*z_esigw3(ji,jj)
551               e3uw_0(ji,jj)=hbatu(ji,jj)*z_esigwu3(ji,jj)
552               e3vw_0(ji,jj)=hbatv(ji,jj)*z_esigwv3(ji,jj)
553            ENDDO
554         ENDDO
555         ! Keep some arrays for next level
556         z_gsigt3m1 = z_gsigt3
557         z_gsi3w3m1 = z_gsi3w3
558 
559         where (e3t_0   (:,:).eq.0.0)  e3t_0(:,:) = 1.0
560         where (e3u_0   (:,:).eq.0.0)  e3u_0(:,:) = 1.0
561         where (e3v_0   (:,:).eq.0.0)  e3v_0(:,:) = 1.0
562         where (e3f_0   (:,:).eq.0.0)  e3f_0(:,:) = 1.0
563         where (e3w_0   (:,:).eq.0.0)  e3w_0(:,:) = 1.0
564         where (e3uw_0  (:,:).eq.0.0)  e3uw_0(:,:) = 1.0
565         where (e3vw_0  (:,:).eq.0.0)  e3vw_0(:,:) = 1.0
566         
567         CALL write_netcdf_3d_vars(jk)
568
569         DO jj = 1, jpj
570            DO ji = 1, jpi
571                  IF( scobot(ji,jj) >= gdept_0(ji,jj) )   mbathy(ji,jj) = MAX( 2, jk )
572                  IF( scobot(ji,jj) == 0.             )   mbathy(ji,jj) = 0
573            END DO
574         END DO
575
576      ENDDO ! End of loop over jk
577
578      DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3, z_gsigw3p1, z_gsigt3m1, z_gsi3w3m1)
579      DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
580
581   END SUBROUTINE s_sf12
582
583   SUBROUTINE s_tanh()
584
585      !!----------------------------------------------------------------------
586      !!                  ***  ROUTINE s_tanh***
587      !!                     
588      !! ** Purpose :   stretch the s-coordinate system
589      !!
590      !! ** Method  :   s-coordinate stretch
591      !!
592      !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
593      !!----------------------------------------------------------------------
594
595      USE utils
596      REAL(wp) ::   zcoeft, zcoefw   ! temporary scalars
597
598      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
599      REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
600
601      ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) )
602      ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
603
604      z_gsigw  = 0.   ;   z_gsigt  = 0.   ;   z_gsi3w  = 0.
605      z_esigt  = 0.   ;   z_esigw  = 0. 
606
607      DO jk = 1, jpk
608        z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5 )
609        z_gsigt(jk) = -fssig( REAL(jk,wp)        )
610      END DO
611      WRITE(*,*) 'z_gsigw 1 jpk    ', z_gsigw(1), z_gsigw(jpk)
612      !
613      ! Coefficients for vertical scale factors at w-, t- levels
614!!gm bug :  define it from analytical function, not like juste bellow....
615!!gm        or betteroffer the 2 possibilities....
616      DO jk = 1, jpk-1
617         z_esigt(jk  ) = z_gsigw(jk+1) - z_gsigw(jk)
618         z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
619      END DO
620      z_esigw( 1 ) = 2. * ( z_gsigt(1  ) - z_gsigw(1  ) ) 
621      z_esigt(jpk) = 2. * ( z_gsigt(jpk) - z_gsigw(jpk) )
622      !
623      ! Coefficients for vertical depth as the sum of e3w scale factors
624      z_gsi3w(1) = 0.5 * z_esigw(1)
625      DO jk = 2, jpk
626         z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
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_3d_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
Note: See TracBrowser for help on using the repository browser.