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.
utils.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/utils.F90 @ 5994

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

Added rx1 to output file as requested by reviewer
File size: 16.4 KB
Line 
1MODULE utils
2
3   USE netcdf
4
5   IMPLICIT NONE
6   PUBLIC             ! allows the acces to par_oce when dom_oce is used
7   !                  ! exception to coding rules... to be suppressed ???
8
9!   PUBLIC dom_oce_alloc
10
11   INTEGER, PARAMETER   :: dp=8 , sp=4, wp=dp
12
13   !! All coordinates
14   !! ---------------
15   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdep3w_0           !: depth of t-points (sum of e3w) (m)
16   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gdept_0, gdepw_0   !: analytical (time invariant) depth at t-w  points (m)
17   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   gphit              !: latitude at t points
18   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3v_0  , e3f_0     !: analytical (time invariant) vertical scale factors at  v-f
19   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3t_0  , e3u_0     !:                                      t-u  points (m)
20   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3vw_0             !: analytical (time invariant) vertical scale factors at  vw
21   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e3w_0  , e3uw_0    !:                                      w-uw points (m)
22
23   !! s-coordinate and hybrid z-s-coordinate
24   !! =----------------======---------------
25   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatv , hbatf      !: ocean depth at the vertical of  v--f
26   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hbatt , hbatu      !:                                 t--u points (m)
27   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   scosrf, scobot     !: ocean surface and bottom topographies
28   !                                                                           !  (if deviating from coordinate surfaces in HYBRID)
29   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hifv  , hiff       !: interface depth between stretching at v--f
30   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hift  , hifu       !: and quasi-uniform spacing             t--u points (m)
31   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rx1                !: Maximum grid stiffness ratio
32
33   !!----------------------------------------------------------------------
34   !! masks, bathymetry
35   !! ---------------------------------------------------------------------
36   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   mbathy             !: number of ocean level (=0, 1, ... , jpk-1)
37   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bathy              !: ocean depth (meters - read from file)
38
39   !! Other variables needed by scoord_gen
40   INTEGER  ::   jpi, jpj, jpk            ! Size of the domain - read from bathy or namelist?
41   INTEGER  ::   ji, jj, jk, jl           ! dummy loop argument
42   INTEGER  ::   iip1, ijp1, iim1, ijm1   ! temporary integers
43   INTEGER  ::   ios                      ! Local integer output status for namelist read and allocation
44   INTEGER,PARAMETER  ::   numnam=8       ! File handle for namelist
45   REAL(wp) ::   zrmax, ztaper   ! temporary scalars
46   REAL(wp) ::   zrfact
47   !
48   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
49   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zenv, ztmp, zmsk, zri, zrj, zhbat
50
51   !Namelist variables
52   REAL(wp) :: rn_jpk, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax, rn_theta
53   REAL(wp) :: rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
54   LOGICAL :: ln_s_sh94, ln_s_sf12, ln_sigcrit, ln_eq_taper
55   CHARACTER(len=50) :: cn_coord_hgr
56
57   NAMELIST/namzgr_sco/rn_jpk, ln_s_sh94, ln_s_sf12, ln_sigcrit, ln_eq_taper, & 
58                  &      cn_coord_hgr, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
59                  &      rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
60
61  ! IDs for output netcdf file
62  INTEGER :: id_x, id_y, id_z
63  INTEGER :: ncout
64  INTEGER, DIMENSION(20) :: var_ids  !Array to contain all variable IDs
65
66   CONTAINS
67
68   INTEGER FUNCTION dom_oce_alloc()
69      !!----------------------------------------------------------------------
70      INTEGER, DIMENSION(4) :: ierr
71      !!----------------------------------------------------------------------
72      ierr(:) = 0
73      !
74      ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), &
75         &      zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj), STAT=ierr(1) )
76         !
77      ALLOCATE( gdep3w_0(jpi,jpj) , e3v_0(jpi,jpj) , e3f_0(jpi,jpj) ,                         &
78         &      gdept_0(jpi,jpj) , e3t_0(jpi,jpj) , e3u_0 (jpi,jpj) ,                         &
79         &      gdepw_0(jpi,jpj) , e3w_0(jpi,jpj) , e3vw_0(jpi,jpj) ,                         &
80         &      gphit(jpi,jpj)   , e3uw_0(jpi,jpj) , STAT=ierr(2) )
81         !
82         !
83      ALLOCATE( hbatv (jpi,jpj) , hbatf (jpi,jpj) ,     &
84         &      hbatt (jpi,jpj) , hbatu (jpi,jpj) ,     &
85         &      scosrf(jpi,jpj) , scobot(jpi,jpj) ,     &
86         &      hifv  (jpi,jpj) , hiff  (jpi,jpj) ,     &
87         &      hift  (jpi,jpj) , hifu  (jpi,jpj) , rx1 (jpi,jpj) , STAT=ierr(3) )
88
89      ALLOCATE( mbathy(jpi,jpj) , STAT=ierr(4) )
90     !
91      dom_oce_alloc = MAXVAL(ierr)
92      !
93   END FUNCTION dom_oce_alloc
94 
95
96   SUBROUTINE read_bathy()
97     !! Read bathymetry from input netcdf file
98     INTEGER :: var_id, ncin
99
100     CALL check_nf90( nf90_open('bathy_meter.nc', NF90_NOWRITE, ncin), 'Error opening bathy_meter.nc file' )
101
102     ! Find the size of the input bathymetry
103     CALL dimlen(ncin, 'lon', jpi)   
104     CALL dimlen(ncin, 'lat', jpj)   
105     
106     ALLOCATE( bathy(jpi, jpj) )
107     
108     ! Read the bathymetry variable from file
109     CALL check_nf90( nf90_inq_varid( ncin, 'Bathymetry', var_id ), 'Cannot get variable ID for bathymetry')
110     CALL check_nf90( nf90_get_var( ncin, var_id, bathy, (/ 1,1 /), (/ jpi, jpj /) ) )
111
112     CALL check_nf90( nf90_close(ncin), 'Error closing bathy.nc file' )
113
114   END SUBROUTINE read_bathy
115
116   SUBROUTINE read_gphit()
117   !! Read gphit from horizontal coordinate file if required
118     INTEGER :: var_id, ncin
119
120     CALL check_nf90( nf90_open(cn_coord_hgr, NF90_NOWRITE, ncin), 'Error opening horizontal coordinate file' )
121
122     ! Read gphit variable from file
123     CALL check_nf90( nf90_inq_varid( ncin, 'gphit', var_id ), 'Cannot get variable ID for bathymetry')
124     CALL check_nf90( nf90_get_var( ncin, var_id, gphit, (/ 1,1 /), (/ jpi, jpj /) ) )
125
126     CALL check_nf90( nf90_close(ncin), 'Error closing horizontal coordinate file' )
127
128   END SUBROUTINE read_gphit
129
130   SUBROUTINE dimlen( ncid, dimname, len )
131     ! Determine the length of dimension dimname
132     INTEGER, INTENT(in)          :: ncid
133     CHARACTER(LEN=*), INTENT(in) :: dimname
134     INTEGER, INTENT(out)         :: len
135     ! Local variables
136     INTEGER :: id_var
137
138     id_var = 1
139     CALL check_nf90( nf90_inq_dimid(ncid, dimname, id_var), 'Dimension not found in file')
140     CALL check_nf90( nf90_inquire_dimension(ncid,id_var,len=len))
141
142   END SUBROUTINE dimlen
143 
144 
145   SUBROUTINE make_coord_file()
146     ! Create new coordinates file and define dimensions and variables ready for
147     ! writing
148     
149
150     !Create the file
151     CALL check_nf90( nf90_create('coord_zgr.nc', NF90_NETCDF4, ncout), 'Could not create output file')
152     !
153     !Define dimensions
154     CALL check_nf90( nf90_def_dim(ncout, 'x', jpi, id_x) )
155     CALL check_nf90( nf90_def_dim(ncout, 'y', jpj, id_y) )
156     CALL check_nf90( nf90_def_dim(ncout, 'z', jpk, id_z) )
157     !
158     !Define variables - include all varibles that would be put into the mesh
159     !mask file
160     CALL check_nf90( nf90_def_var(ncout, 'gdept_0', nf90_double, (/id_x, id_y,id_z/), var_ids(1)) )
161     CALL check_nf90( nf90_def_var(ncout, 'gdepw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(2)) )
162     CALL check_nf90( nf90_def_var(ncout, 'gdep3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(3)) )
163     CALL check_nf90( nf90_def_var(ncout, 'e3f_0', nf90_double, (/id_x, id_y,id_z/), var_ids(4)) )
164     CALL check_nf90( nf90_def_var(ncout, 'e3t_0', nf90_double, (/id_x, id_y,id_z/), var_ids(5)) )
165     CALL check_nf90( nf90_def_var(ncout, 'e3u_0', nf90_double, (/id_x, id_y,id_z/), var_ids(6)) )
166     CALL check_nf90( nf90_def_var(ncout, 'e3v_0', nf90_double, (/id_x, id_y,id_z/), var_ids(7)) )
167     CALL check_nf90( nf90_def_var(ncout, 'e3w_0', nf90_double, (/id_x, id_y,id_z/), var_ids(8)) )
168     CALL check_nf90( nf90_def_var(ncout, 'e3uw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(9)) )
169     CALL check_nf90( nf90_def_var(ncout, 'e3vw_0', nf90_double, (/id_x, id_y,id_z/), var_ids(10)) )
170     ! 2D fields
171     CALL check_nf90( nf90_def_var(ncout, 'mbathy', nf90_double, (/id_x, id_y/), var_ids(11)) )
172     CALL check_nf90( nf90_def_var(ncout, 'hbatt', nf90_double, (/id_x, id_y/), var_ids(12)) )
173     CALL check_nf90( nf90_def_var(ncout, 'hbatu', nf90_double, (/id_x, id_y/), var_ids(13)) )
174     CALL check_nf90( nf90_def_var(ncout, 'hbatv', nf90_double, (/id_x, id_y/), var_ids(14)) )
175     CALL check_nf90( nf90_def_var(ncout, 'hbatf', nf90_double, (/id_x, id_y/), var_ids(15)) )
176     CALL check_nf90( nf90_def_var(ncout, 'rx1', nf90_double, (/id_x, id_y/), var_ids(16)) )
177
178     
179     ! End define mode
180     CALL check_nf90( nf90_enddef(ncout) )
181     
182     WRITE(*,*) 'Opened coord_zgr.nc file and defined variables'
183
184   END SUBROUTINE make_coord_file
185
186   SUBROUTINE write_netcdf_2d_vars()
187
188     CALL check_nf90( nf90_put_var(ncout, var_ids(11), mbathy, (/ 1,1 /), (/ jpi, jpj /) ) )
189     CALL check_nf90( nf90_put_var(ncout, var_ids(12), hbatt, (/ 1,1 /), (/ jpi, jpj /) ) )
190     CALL check_nf90( nf90_put_var(ncout, var_ids(13), hbatu, (/ 1,1 /), (/ jpi, jpj /) ) )
191     CALL check_nf90( nf90_put_var(ncout, var_ids(14), hbatv, (/ 1,1 /), (/ jpi, jpj /) ) )
192     CALL check_nf90( nf90_put_var(ncout, var_ids(15), hbatf, (/ 1,1 /), (/ jpi, jpj /) ) )
193     CALL check_nf90( nf90_put_var(ncout, var_ids(16), rx1, (/ 1,1 /), (/ jpi, jpj /) ) )
194
195   END SUBROUTINE write_netcdf_2d_vars
196
197   SUBROUTINE write_netcdf_3d_vars(kk)
198   ! Write  variables to the netcdf file at level kk
199     INTEGER, INTENT(in) :: kk
200
201     CALL check_nf90( nf90_put_var(ncout, var_ids(1), gdept_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
202     CALL check_nf90( nf90_put_var(ncout, var_ids(2), gdepw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
203     CALL check_nf90( nf90_put_var(ncout, var_ids(3), gdep3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
204     CALL check_nf90( nf90_put_var(ncout, var_ids(4), e3f_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
205     CALL check_nf90( nf90_put_var(ncout, var_ids(5), e3t_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
206     CALL check_nf90( nf90_put_var(ncout, var_ids(6), e3u_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
207     CALL check_nf90( nf90_put_var(ncout, var_ids(7), e3v_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
208     CALL check_nf90( nf90_put_var(ncout, var_ids(8), e3w_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
209     CALL check_nf90( nf90_put_var(ncout, var_ids(9), e3uw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
210     CALL check_nf90( nf90_put_var(ncout, var_ids(10), e3vw_0, (/ 1,1,kk /), (/ jpi, jpj,1 /) ) )
211
212   END SUBROUTINE write_netcdf_3d_vars
213
214   SUBROUTINE check_nf90( istat, message )
215      !Check for netcdf errors
216      INTEGER, INTENT(in) :: istat
217      CHARACTER(LEN=*), INTENT(in), OPTIONAL :: message
218
219      IF (istat /= nf90_noerr) THEN
220         WRITE(*,*) 'ERROR! : '//TRIM(nf90_strerror(istat))
221         IF ( PRESENT(message) ) THEN ; WRITE(*,*) message ; ENDIF
222         STOP
223      ENDIF
224
225   END SUBROUTINE check_nf90
226   FUNCTION fssig( pk ) RESULT( pf )
227      !!----------------------------------------------------------------------
228      !!                 ***  ROUTINE fssig ***
229      !!       
230      !! ** Purpose :   provide the analytical function in s-coordinate
231      !!         
232      !! ** Method  :   the function provide the non-dimensional position of
233      !!                T and W (i.e. between 0 and 1)
234      !!                T-points at integer values (between 1 and jpk)
235      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
236      !!----------------------------------------------------------------------
237!     USE utils, ONLY : wp,rn_theta,rn_thetb,jpk
238      IMPLICIT NONE
239      REAL(wp), INTENT(in) ::   pk   ! continuous "k" coordinate
240      REAL(wp)             ::   pf   ! sigma value
241      !!----------------------------------------------------------------------
242      !
243      pf =   (   TANH( rn_theta * ( -(pk-0.5) / REAL(jpk-1) + rn_thetb )  )   &
244         &     - TANH( rn_thetb * rn_theta                                )  )   &
245         & * (   COSH( rn_theta                           )                      &
246         &     + COSH( rn_theta * ( 2. * rn_thetb - 1. ) )  )              &
247         & / ( 2. * SINH( rn_theta ) )
248      !
249   END FUNCTION fssig
250
251   FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
252      !!----------------------------------------------------------------------
253      !!                 ***  ROUTINE fssig1 ***
254      !!
255      !! ** Purpose :   provide the Song and Haidvogel version of the analytical function in s-coordinate
256      !!
257      !! ** Method  :   the function provides the non-dimensional position of
258      !!                T and W (i.e. between 0 and 1)
259      !!                T-points at integer values (between 1 and jpk)
260      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
261      !!----------------------------------------------------------------------
262!     USE utils, ONLY : wp, jpk, rn_theta
263      IMPLICIT NONE
264      REAL(wp), INTENT(in) ::   pk1   ! continuous "k" coordinate
265      REAL(wp), INTENT(in) ::   pbb   ! Stretching coefficient
266      REAL(wp)             ::   pf1   ! sigma value
267      !!----------------------------------------------------------------------
268      !
269      IF ( rn_theta == 0 ) then      ! uniform sigma
270         pf1 = - ( pk1 - 0.5 ) / REAL( jpk-1 )
271      ELSE                        ! stretched sigma
272         pf1 =   ( 1. - pbb ) * ( SINH( rn_theta*(-(pk1-0.5)/REAL(jpk-1)) ) ) / SINH( rn_theta )              &
273            &  + pbb * (  (TANH( rn_theta*( (-(pk1-0.5)/REAL(jpk-1)) + 0.5) ) - TANH( 0.5 * rn_theta )  )  &
274            &        / ( 2. * TANH( 0.5 * rn_theta ) )  )
275      ENDIF
276      !
277   END FUNCTION fssig1
278
279   FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
280      !!----------------------------------------------------------------------
281      !!                 ***  ROUTINE fgamma  ***
282      !!
283      !! ** Purpose :   provide analytical function for the s-coordinate
284      !!
285      !! ** Method  :   the function provides the non-dimensional position of
286      !!                T and W (i.e. between 0 and 1)
287      !!                T-points at integer values (between 1 and jpk)
288      !!                W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
289      !!
290      !!                This method allows the maintenance of fixed surface and or
291      !!                bottom cell resolutions (cf. geopotential coordinates)
292      !!                within an analytically derived stretched S-coordinate framework.
293      !!
294      !! Reference  :   Siddorn and Furner, in prep
295      !!----------------------------------------------------------------------
296!     USE utils, ONLY : jpk,wp,rn_alpha
297      IMPLICIT NONE
298      REAL(wp), INTENT(in   ) ::   pk1           ! continuous "k" coordinate
299      REAL(wp)                ::   p_gamma       ! stretched coordinate
300      REAL(wp), INTENT(in   ) ::   pzb           ! Bottom box depth
301      REAL(wp), INTENT(in   ) ::   pzs           ! surface box depth
302      REAL(wp), INTENT(in   ) ::   psmth         ! Smoothing parameter
303      REAL(wp)                ::   za1,za2,za3   ! local variables
304      REAL(wp)                ::   zn1,zn2       ! local variables
305      REAL(wp)                ::   za,zb,zx      ! local variables
306      !!----------------------------------------------------------------------
307      !
308
309      zn1  =  1./(jpk-1.)
310      zn2  =  1. -  zn1
311
312      za1 = (rn_alpha+2.0)*zn1**(rn_alpha+1.0)-(rn_alpha+1.0)*zn1**(rn_alpha+2.0) 
313      za2 = (rn_alpha+2.0)*zn2**(rn_alpha+1.0)-(rn_alpha+1.0)*zn2**(rn_alpha+2.0)
314      za3 = (zn2**3.0 - za2)/( zn1**3.0 - za1)
315     
316      za = pzb - za3*(pzs-za1)-za2
317      za = za/( zn2-0.5*(za2+zn2**2.0) - za3*(zn1-0.5*(za1+zn1**2.0) ) )
318      zb = (pzs - za1 - za*( zn1-0.5*(za1+zn1**2.0 ) ) ) / (zn1**3.0 - za1)
319      zx = 1.0-za/2.0-zb
320
321      p_gamma = za*(pk1*(1.0-pk1/2.0))+zb*pk1**3.0 +  &
322                  & zx*( (rn_alpha+2.0)*pk1**(rn_alpha+1.0)- &
323                  &      (rn_alpha+1.0)*pk1**(rn_alpha+2.0) )
324      p_gamma = p_gamma*psmth+pk1*(1.0-psmth)
325
326      !
327   END FUNCTION fgamma
328
329
330END MODULE utils
Note: See TracBrowser for help on using the repository browser.