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.
vgrid.f90 in NEMO/trunk/tools/SIREN/src – NEMO

source: NEMO/trunk/tools/SIREN/src/vgrid.f90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 6 years ago

Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure
Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing
Standardisation of routines header with version 4.0 and year 2018
Fix for some broken symlinks

File size: 32.0 KB
Line 
1!----------------------------------------------------------------------
2! NEMO system team, System and Interface for oceanic RElocable Nesting
3!----------------------------------------------------------------------
4!
5! MODULE: vgrid
6!
7! DESCRIPTION:
8!> @brief This module manage vertical grid.
9!>
10!> @details
11!>    to set the depth of model levels and the resulting vertical scale
12!> factors:<br/>
13!> @code
14!>    CALL vgrid_zgr_z(dd_gdepw(:), dd_gdept(:), dd_e3w(:), dd_e3t(:),
15!>                     dd_ppkth, dd_ppkth2, dd_ppacr, dd_ppacr2,
16!>                     dd_ppdzmin, dd_pphmax,
17!>                     dd_ppa0, dd_ppa1, dd_ppa2, dd_ppsur)
18!> @endcode
19!>       - dd_gdepw is array of depth value on W point
20!>       - dd_gdept is array of depth value on T point
21!>       - dd_e3w   is array of vertical mesh size on W point
22!>       - dd_e3t   is array of vertical mesh size on T point
23!>       - dd_ppkth              see NEMO documentation
24!>       - dd_ppkth2             see NEMO documentation
25!>       - dd_ppacr              see NEMO documentation
26!>       - dd_ppdzmin            see NEMO documentation
27!>       - dd_pphmax             see NEMO documentation
28!>       - dd_ppa1               see NEMO documentation
29!>       - dd_ppa2               see NEMO documentation
30!>       - dd_ppa0               see NEMO documentation
31!>       - dd_ppsur              see NEMO documentation
32!>   
33!>
34!>    to set the depth and vertical scale factor in partial step z-coordinate
35!>  case:<br/>
36!> @code
37!>    CALL vgrid_zgr_zps(id_mbathy(:,:), dd_bathy(:,:), id_jpkmax, dd_gdepw(:),
38!>                       dd_e3t(:), dd_e3zps_min, dd_e3zps_rat)
39!> @endcode
40!>       - id_mbathy is array of bathymetry level
41!>       - dd_bathy  is array of bathymetry
42!>       - id_jpkmax is the maximum number of level to be used
43!>       - dd_gdepw  is array of vertical mesh size on W point
44!>       - dd_e3t    is array of vertical mesh size on T point
45!>       - dd_e3zps_min    see NEMO documentation
46!>       - dd_e3zps_rat    see NEMO documentation
47!>
48!>    to check the bathymetry in levels:<br/>
49!> @code
50!>    CALL vgrid_zgr_bat_ctl(id_mbathy, id_jpkmax, id_jpk)
51!> @endcode
52!>       - id_mbathy is array of bathymetry level
53!>       - id_jpkmax is the maximum number of level to be used
54!>       - id_jpk    is the number of level
55!>   
56!>    to compute bathy level in T,U,V,F point from  Bathymetry file:<br/>
57!> @code
58!>    tl_level(:)=vgrid_get_level(td_bathy, [cd_namelist,] [td_dom,] [id_nlevel])
59!> @endcode
60!>       - td_bathy is Bathymetry file structure
61!>       - cd_namelist is namelist [optional]
62!>       - td_dom is domain structure [optional]
63!>       - id_nlevel is number of lelvel to be used [optional]
64!>   
65!> @author
66!> J.Paul
67! REVISION HISTORY:
68!> @date November, 2013 - Initial Version
69!> @date Spetember, 2014
70!> - add header
71!> @date June, 2015 - update subroutine with NEMO 3.6
72!>
73!> @note Software governed by the CeCILL licence     (./LICENSE)
74!----------------------------------------------------------------------
75MODULE vgrid
76   USE netcdf                          ! nf90 library
77   USE kind                            ! F90 kind parameter
78   USE fct                             ! basic usefull function
79   USE global                          ! global parameter
80   USE phycst                          ! physical constant
81   USE logger                          ! log file manager
82   USE file                            ! file manager
83   USE var                             ! variable manager
84   USE dim                             ! dimension manager
85   USE dom                             ! domain manager
86   USE grid                            ! grid manager
87   USE iom                             ! I/O manager
88   USE mpp                             ! MPP manager
89   USE iom_mpp                         ! I/O MPP manager
90   IMPLICIT NONE
91   ! NOTE_avoid_public_variables_if_possible
92
93   ! type and variable
94
95   ! function and subroutine
96   PUBLIC :: vgrid_zgr_z 
97   PUBLIC :: vgrid_zgr_zps
98   PUBLIC :: vgrid_zgr_bat_ctl
99   PUBLIC :: vgrid_get_level
100
101CONTAINS
102   !-------------------------------------------------------------------
103   !> @brief This subroutine set the depth of model levels and the resulting
104   !>      vertical scale factors.
105   !
106   !> @details
107   !> ** Method  :   z-coordinate system (use in all type of coordinate)
108   !>        The depth of model levels is defined from an analytical
109   !>      function the derivative of which gives the scale factors.
110   !>        both depth and scale factors only depend on k (1d arrays). <\br>
111   !>              w-level: gdepw  = fsdep(k)                           <\br>
112   !>                       e3w(k) = dk(fsdep)(k)     = fse3(k)         <\br>
113   !>              t-level: gdept  = fsdep(k+0.5)                       <\br>
114   !>                       e3t(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)     <\br>
115   !>
116   !> ** Action  : - gdept, gdepw : depth of T- and W-point (m)          <\br>
117   !>              -  e3t, e3w    : scale factors at T- and W-levels (m) <\br>
118   !>
119   !> @author G. Madec
120   !> @date Marsh,2008 - F90: Free form and module
121   !
122   !> @note Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
123   !>
124   !> @param[inout] dd_gdepw
125   !> @param[inout] dd_gedpt
126   !> @param[inout] dd_e3w
127   !> @param[inout] dd_e2t
128   !> @param[in] dd_ppkth
129   !> @param[in] dd_ppkth2
130   !> @param[in] dd_ppacr
131   !> @param[in] dd_ppacr2
132   !> @param[in] dd_ppdzmin
133   !> @param[in] dd_pphmax
134   !> @param[in] dd_ppa1
135   !> @param[in] dd_ppa2
136   !> @param[in] dd_ppa0
137   !> @param[in] dd_ppsur
138   !-------------------------------------------------------------------
139   SUBROUTINE vgrid_zgr_z( dd_gdepw, dd_gdept, dd_e3w, dd_e3t,          &
140   &                       dd_e3w_1d, dd_e3t_1d,                        &
141   &                       dd_ppkth, dd_ppkth2, dd_ppacr, dd_ppacr2,    &
142   &                       dd_ppdzmin, dd_pphmax,                       &
143   &                       dd_ppa0, dd_ppa1, dd_ppa2, dd_ppsur )
144      IMPLICIT NONE
145      ! Argument     
146      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_gdepw
147      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_gdept
148      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3w
149      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3t
150      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3w_1d
151      REAL(dp), DIMENSION(:), INTENT(INOUT) :: dd_e3t_1d
152
153      REAL(dp)              , INTENT(IN   ) :: dd_ppkth
154      REAL(dp)              , INTENT(IN   ) :: dd_ppkth2
155      REAL(dp)              , INTENT(IN   ) :: dd_ppacr
156      REAL(dp)              , INTENT(IN   ) :: dd_ppacr2
157
158      REAL(dp)              , INTENT(IN   ) :: dd_ppdzmin
159      REAL(dp)              , INTENT(IN   ) :: dd_pphmax
160      REAL(dp), PARAMETER                   :: dp_pp_to_be_computed = NF90_FILL_DOUBLE
161
162      REAL(dp)              , INTENT(IN   ) :: dd_ppa0
163      REAL(dp)              , INTENT(IN   ) :: dd_ppa1
164      REAL(dp)              , INTENT(IN   ) :: dd_ppa2
165      REAL(dp)              , INTENT(IN   ) :: dd_ppsur
166
167      ! local variable
168      REAL(dp) :: dl_zkth
169      REAL(dp) :: dl_zkth2
170      REAL(dp) :: dl_zdzmin
171      REAL(dp) :: dl_zhmax
172      REAL(dp) :: dl_zacr
173      REAL(dp) :: dl_zacr2
174
175      REAL(dp) :: dl_ppacr
176      REAL(dp) :: dl_ppacr2
177
178      REAL(dp) :: dl_za0
179      REAL(dp) :: dl_za1
180      REAL(dp) :: dl_za2
181      REAL(dp) :: dl_zsur
182      REAL(dp) :: dl_zw
183      REAL(dp) :: dl_zt
184
185      INTEGER(i4) :: il_jpk
186
187      ! loop indices
188      INTEGER(i4) :: jk
189      !----------------------------------------------------------------
190
191       dl_ppacr = dd_ppacr
192       IF( dd_ppa1 == 0._dp ) dl_ppacr =1.0
193       dl_ppacr2= dd_ppacr2
194       IF( dd_ppa2 == 0._dp ) dl_ppacr2=1.0
195
196      ! Set variables from parameters
197      ! ------------------------------
198       dl_zkth   = dd_ppkth       ;   dl_zacr  = dl_ppacr
199       dl_zdzmin = dd_ppdzmin     ;   dl_zhmax = dd_pphmax
200       dl_zkth2  = dd_ppkth2      ;   dl_zacr2 = dl_ppacr2
201
202       il_jpk = SIZE(dd_gdepw(:))
203
204      ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
205      !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
206      !
207       IF(  dd_ppa1  == dp_pp_to_be_computed  .AND.  &
208         &  dd_ppa0  == dp_pp_to_be_computed  .AND.  &
209         &  dd_ppsur == dp_pp_to_be_computed           ) THEN
210         dl_za1 = ( dl_zdzmin - dl_zhmax / REAL((il_jpk-1),dp) ) &
211             &     / ( TANH((1-dl_zkth)/dl_zacr) - dl_zacr/REAL((il_jpk-1),dp) &
212             &     * (  LOG( COSH( (REAL(il_jpk,dp) - dl_zkth) / dl_zacr) )    &
213             &        - LOG( COSH(( 1.0  - dl_zkth) / dl_zacr) ) )  )
214
215         dl_za0  = dl_zdzmin - dl_za1 * TANH( (1.0-dl_zkth) / dl_zacr )
216         dl_zsur = - dl_za0 - dl_za1 * dl_zacr * LOG( COSH( (1-dl_zkth) / dl_zacr )  )
217
218       ELSE
219         dl_za1 = dd_ppa1 ;       dl_za0 = dd_ppa0 ;  dl_zsur = dd_ppsur
220         dl_za2 = dd_ppa2
221       ENDIF
222
223      ! Reference z-coordinate (depth - scale factor at T- and W-points)
224      ! ======================
225      IF(  dd_ppkth == 0. )THEN            !  uniform vertical grid       
226
227         dl_za1 = dl_zhmax/REAL((il_jpk-1),dp) 
228         DO jk = 1, il_jpk
229            dl_zw = REAL(jk,dp)
230            dl_zt = REAL(jk,dp) + 0.5_dp
231            dd_gdepw(jk) = ( dl_zw - 1.0 ) * dl_za1
232            dd_gdept(jk) = ( dl_zt - 1.0 ) * dl_za1
233            dd_e3w  (jk) =  dl_za1
234            dd_e3t  (jk) =  dl_za1
235         END DO
236
237      ELSE
238
239         DO jk = 1, il_jpk
240            dl_zw = REAL( jk,dp)
241            dl_zt = REAL( jk,dp) + 0.5_dp
242            dd_gdepw(jk) = ( dl_zsur + dl_za0 * dl_zw + &
243            &                dl_za1 * dl_zacr * LOG( COSH( (dl_zw-dl_zkth)/dl_zacr ) ) + &
244            &                dl_za2 * dl_zacr2* LOG( COSH( (dl_zw-dl_zkth2)/dl_zacr2 ) ) )
245            dd_gdept(jk) = ( dl_zsur + dl_za0 * dl_zt + &
246            &                dl_za1 * dl_zacr * LOG( COSH( (dl_zt-dl_zkth)/dl_zacr ) ) + &
247            &                dl_za2 * dl_zacr2* LOG( COSH( (dl_zt-dl_zkth2)/dl_zacr2 ) ) )
248            dd_e3w  (jk) =             dl_za0 + &
249            &                          dl_za1 * TANH(      (dl_zw-dl_zkth)/dl_zacr   ) + &
250            &                          dl_za2 * TANH(      (dl_zw-dl_zkth2)/dl_zacr2 )
251            dd_e3t  (jk) =             dl_za0 + &
252            &                          dl_za1 * TANH(      (dl_zt-dl_zkth)/dl_zacr   ) + &
253            &                          dl_za2 * TANH(      (dl_zt-dl_zkth2)/dl_zacr2 )
254         END DO
255         dd_gdepw(1) = 0.e0   ! force first w-level to be exactly at zero
256
257      ENDIF
258
259   ! need to be like this to compute the pressure gradient with ISF.
260   ! If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
261   ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
262      DO jk = 1, il_jpk-1
263         dd_e3t_1d(jk) = dd_gdepw(jk+1)-dd_gdepw(jk) 
264      END DO
265      dd_e3t_1d(il_jpk) = dd_e3t_1d(il_jpk-1) ! we don't care because this level is masked in NEMO
266
267      DO jk = 2, il_jpk
268         dd_e3w_1d(jk) = dd_gdept(jk) - dd_gdept(jk-1) 
269      END DO
270      dd_e3w_1d(1  ) = 2._dp * (dd_gdept(1) - dd_gdepw(1))
271
272      ! Control and  print
273      ! ==================
274
275      DO jk = 1, il_jpk
276         IF( dd_e3w(jk)  <= 0. .OR. dd_e3t(jk)  <= 0. )then
277            CALL logger_debug("VGRID ZGR Z: e3w or e3t <= 0 ")
278         ENDIF   
279
280         IF( dd_e3w_1d(jk)  <= 0. .OR. dd_e3t_1d(jk)  <= 0. )then
281            CALL logger_debug("VGRID ZGR Z: e3w_1d or e3t_1d <= 0 ")
282         ENDIF
283
284         IF( dd_gdepw(jk) < 0. .OR. dd_gdept(jk) < 0. )then
285            CALL logger_debug("VGRID ZGR Z: gdepw or gdept < 0 ")
286         ENDIF
287      END DO
288
289   END SUBROUTINE vgrid_zgr_z
290   !-------------------------------------------------------------------
291   !> @brief This subroutine
292   !>
293   !> @todo add subroutine description
294   !>
295   !> @param[inout] dd_bathy
296   !> @param[in] dd_gdepw
297   !> @param[in] dd_hmin
298   !> @param[in] dd_fill
299   !-------------------------------------------------------------------
300   SUBROUTINE vgrid_zgr_bat( dd_bathy, dd_gdepw, dd_hmin, dd_fill )
301      IMPLICIT NONE
302      ! Argument
303      REAL(dp), DIMENSION(:,:), INTENT(INOUT) :: dd_bathy 
304      REAL(dp), DIMENSION(:)  , INTENT(IN   ) :: dd_gdepw 
305      REAL(dp)                , INTENT(IN   ) :: dd_hmin
306      REAL(dp)                , INTENT(IN   ), OPTIONAL :: dd_fill
307
308      ! local
309      INTEGER(i4) :: il_jpk
310     
311      REAL(dp)    :: dl_hmin
312      REAL(dp)    :: dl_fill
313
314      ! loop indices
315      INTEGER(i4) :: jk
316      !----------------------------------------------------------------
317      il_jpk = SIZE(dd_gdepw(:))
318
319      dl_fill=0._dp
320      IF( PRESENT(dd_fill) ) dl_fill=dd_fill
321
322      IF( dd_hmin < 0._dp ) THEN
323         jk = - INT( dd_hmin )     ! from a nb of level
324      ELSE
325         jk = MINLOC( dd_gdepw, mask = dd_gdepw > dd_hmin, dim = 1 )  ! from a depth
326      ENDIF
327     
328      dl_hmin = dd_gdepw(jk+1) ! minimum depth = ik+1 w-levels
329      WHERE( dd_bathy(:,:) <= 0._wp .OR. dd_bathy(:,:) == dl_fill )
330         dd_bathy(:,:) = dl_fill                         ! min=0     over the lands
331      ELSE WHERE
332         dd_bathy(:,:) = MAX(  dl_hmin , dd_bathy(:,:)  )   ! min=dl_hmin over the oceans
333      END WHERE
334      WRITE(*,*) 'Minimum ocean depth: ', dl_hmin, ' minimum number of ocean levels : ', jk     
335
336   END SUBROUTINE vgrid_zgr_bat
337   !-------------------------------------------------------------------
338   !> @brief This subroutine set the depth and vertical scale factor in partial step
339   !>      z-coordinate case
340   !
341   !> @details
342   !> ** Method  :   Partial steps : computes the 3D vertical scale factors
343   !>      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
344   !>      a partial step representation of bottom topography.
345   !>
346   !>        The reference depth of model levels is defined from an analytical
347   !>      function the derivative of which gives the reference vertical
348   !>      scale factors.
349   !>      From  depth and scale factors reference, we compute there new value
350   !>      with partial steps  on 3d arrays ( i, j, k ).
351   !>
352   !>      w-level:
353   !>          - gdepw_ps(i,j,k)  = fsdep(k)
354   !>          - e3w_ps(i,j,k) = dk(fsdep)(k)     = fse3(i,j,k)
355   !>      t-level:
356   !>          - gdept_ps(i,j,k)  = fsdep(k+0.5)
357   !>          - e3t_ps(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
358   !>
359   !>      With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
360   !>      we find the mbathy index of the depth at each grid point.
361   !>      This leads us to three cases:
362   !>          - bathy = 0 => mbathy = 0
363   !>          - 1 < mbathy < jpkm1   
364   !>          - bathy > gdepw(jpk) => mbathy = jpkm1 
365   !>
366   !>      Then, for each case, we find the new depth at t- and w- levels
367   !>      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
368   !>      and f-points.
369   !>
370   !>        This routine is given as an example, it must be modified
371   !>      following the user s desiderata. nevertheless, the output as
372   !>      well as the way to compute the model levels and scale factors
373   !>      must be respected in order to insure second order accuracy
374   !>      schemes.
375   !>
376   !>  @warning
377   !>         - gdept, gdepw and e3 are positives
378   !>         - gdept_ps, gdepw_ps and e3_ps are positives
379   !>
380   !> @author A. Bozec, G. Madec
381   !> @date February, 2009 - F90: Free form and module
382   !> @date February, 2009
383   !> - A. de Miranda : rigid-lid + islands
384   !>
385   !> @note Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
386   !>
387   !> @param[inout] id_mbathy
388   !> @param[inout] dd_bathy
389   !> @param[inout] id_jpkmax
390   !> @param[in] dd_gdepw
391   !> @param[in] dd_e3t
392   !> @param[in] dd_e3zps_min
393   !> @param[in] dd_e3zps_rat
394   !> @param[in] dd_fill
395   !-------------------------------------------------------------------
396   SUBROUTINE vgrid_zgr_zps( id_mbathy, dd_bathy, id_jpkmax, &
397   &                          dd_gdepw, dd_e3t,               &
398   &                          dd_e3zps_min, dd_e3zps_rat,     &
399   &                          dd_fill )
400      IMPLICIT NONE
401      ! Argument     
402      INTEGER(i4), DIMENSION(:,:), INTENT(  OUT) :: id_mbathy
403      REAL(dp)   , DIMENSION(:,:), INTENT(INOUT) :: dd_bathy
404      INTEGER(i4)                , INTENT(INOUT) :: id_jpkmax
405      REAL(dp)   , DIMENSION(:)  , INTENT(IN   ) :: dd_gdepw
406      REAL(dp)   , DIMENSION(:)  , INTENT(IN   ) :: dd_e3t
407      REAL(dp)                   , INTENT(IN   ) :: dd_e3zps_min
408      REAL(dp)                   , INTENT(IN   ) :: dd_e3zps_rat
409      REAL(dp)                   , INTENT(IN   ), OPTIONAL :: dd_fill
410
411      ! local variable
412      REAL(dp) :: dl_zmax     ! Maximum depth
413      !REAL(dp) :: dl_zmin     ! Minimum depth
414      REAL(dp) :: dl_zdepth   ! Ajusted ocean depth to avoid too small e3t
415      REAL(dp) :: dl_fill     
416
417      INTEGER(i4) :: il_jpk
418      INTEGER(i4) :: il_jpkm1
419      INTEGER(i4) :: il_jpiglo
420      INTEGER(i4) :: il_jpjglo
421
422      ! loop indices
423      INTEGER(i4) :: ji
424      INTEGER(i4) :: jj
425      INTEGER(i4) :: jk
426      !----------------------------------------------------------------
427
428      il_jpk=SIZE(dd_gdepw(:))
429      il_jpiglo=SIZE(id_mbathy(:,:),DIM=1)
430      il_jpjglo=SIZE(id_mbathy(:,:),DIM=2)
431
432      dl_fill=0._dp
433      IF( PRESENT(dd_fill) ) dl_fill=dd_fill
434
435      ! Initialization of constant
436      dl_zmax = dd_gdepw(il_jpk) + dd_e3t(il_jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
437
438      ! bounded value of bathy (min already set at the end of zgr_bat)
439      WHERE( dd_bathy(:,:) /= dl_fill )
440         dd_bathy(:,:) = MIN( dl_zmax ,  dd_bathy(:,:) )
441      END WHERE
442
443      ! bathymetry in level (from bathy_meter)
444      ! ===================
445      il_jpkm1=il_jpk-1
446      ! initialize mbathy to the maximum ocean level available
447      id_mbathy(:,:) = il_jpkm1
448
449      ! storage of land and island's number (zera and negative values) in mbathy
450      DO jj = 1, il_jpjglo
451         DO ji= 1, il_jpiglo
452            IF( dd_bathy(ji,jj) <= 0._dp )THEN
453               id_mbathy(ji,jj) = INT(dd_bathy(ji,jj),i4)
454            ELSEIF( dd_bathy(ji,jj) == dl_fill )THEN
455               id_mbathy(ji,jj) = 0_i4
456            ENDIF
457         END DO
458      END DO
459
460      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
461      ! find the number of ocean levels such that the last level thickness
462      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t (where
463      ! e3t is the reference level thickness
464
465      DO jk = il_jpkm1, 1, -1
466         dl_zdepth = dd_gdepw(jk) + MIN( dd_e3zps_min, dd_e3t(jk)*dd_e3zps_rat )
467
468         DO jj = 1, il_jpjglo
469            DO ji = 1, il_jpiglo
470               IF( dd_bathy(ji,jj) /= dl_fill )THEN
471                  IF( 0. < dd_bathy(ji,jj) .AND. &
472                  &       dd_bathy(ji,jj) <= dl_zdepth ) id_mbathy(ji,jj) = jk-1
473               ENDIF
474            END DO
475         END DO
476      END DO
477
478      ! ================
479      ! Bathymetry check
480      ! ================
481
482      CALL vgrid_zgr_bat_ctl( id_mbathy, id_jpkmax, il_jpk)
483
484   END SUBROUTINE vgrid_zgr_zps
485   !-------------------------------------------------------------------
486   !> @brief This subroutine check the bathymetry in levels
487   !
488   !> @details
489   !> ** Method  :   The array mbathy is checked to verified its consistency
490   !>      with the model options. in particular:
491   !>            mbathy must have at least 1 land grid-points (mbathy<=0)
492   !>                  along closed boundary.
493   !>            mbathy must be cyclic IF jperio=1.
494   !>            mbathy must be lower or equal to jpk-1.
495   !>            isolated ocean grid points are suppressed from mbathy
496   !>                  since they are only connected to remaining
497   !>                  ocean through vertical diffusion.
498   !>      C A U T I O N : mbathy will be modified during the initializa-
499   !>      tion phase to become the number of non-zero w-levels of a water
500   !>      column, with a minimum value of 1.
501   !>
502   !> ** Action  : - update mbathy: level bathymetry (in level index)
503   !>              - update bathy : meter bathymetry (in meters)
504   !>
505   !> @author G.Madec
506   !> @date Marsh, 2008 - Original code
507   !
508   !> @param[in] id_mbathy
509   !> @param[in] id_jpkmax
510   !> @param[in] id_jpk
511   !-------------------------------------------------------------------
512   SUBROUTINE vgrid_zgr_bat_ctl( id_mbathy, id_jpkmax, id_jpk)
513      IMPLICIT NONE
514      ! Argument     
515      INTEGER(i4), DIMENSION(:,:), INTENT(INOUT) :: id_mbathy
516      INTEGER(i4)                , INTENT(INOUT) :: id_jpkmax
517      INTEGER(i4)                , INTENT(INOUT) :: id_jpk
518
519      ! local variable
520      INTEGER(i4) :: il_jpiglo
521      INTEGER(i4) :: il_jpjglo
522
523      INTEGER(i4) :: il_icompt
524      INTEGER(i4) :: il_ibtest
525      INTEGER(i4) :: il_ikmax
526      INTEGER(i4) :: il_jpkm1
527
528      INTEGER(i4) :: il_jim
529      INTEGER(i4) :: il_jip
530      INTEGER(i4) :: il_jjm
531      INTEGER(i4) :: il_jjp
532
533      ! loop indices
534      INTEGER(i4) :: jl
535      INTEGER(i4) :: ji
536      INTEGER(i4) :: jj
537      !----------------------------------------------------------------
538
539      il_jpiglo=SIZE(id_mbathy(:,:),DIM=1)
540      il_jpjglo=SIZE(id_mbathy(:,:),DIM=2)
541
542      ! ================
543      ! Bathymetry check
544      ! ================
545
546      ! suppress isolated ocean grid points'
547      il_icompt = 0
548
549      DO jl = 1, 2
550         DO jj = 1, il_jpjglo
551            DO ji = 1, il_jpiglo
552               il_jim=max(ji-1,1) ; il_jip=min(ji+1,il_jpiglo)
553               il_jjm=max(jj-1,1) ; il_jjp=min(jj+1,il_jpjglo)
554
555               if(il_jim==ji) il_jim=il_jip ; if(il_jip==ji) il_jip=il_jim
556               if(il_jjm==jj) il_jjm=il_jjp ; if(il_jjp==jj) il_jjp=il_jjm
557
558               il_ibtest = MAX( id_mbathy(il_jim,jj), id_mbathy(il_jip,jj),  &
559                                id_mbathy(ji,il_jjm),id_mbathy(ji,il_jjp) )
560
561               IF( il_ibtest < id_mbathy(ji,jj) ) THEN
562                  id_mbathy(ji,jj) = il_ibtest
563                  il_icompt = il_icompt + 1
564               ENDIF
565            END DO
566         END DO
567
568      END DO
569      IF( il_icompt == 0 ) THEN
570         CALL logger_info("VGRID ZGR BAT CTL:  no isolated ocean grid points")
571      ELSE
572         CALL logger_info("VGRID ZGR BAT CTL:"//TRIM(fct_str(il_icompt))//&
573         &              " ocean grid points suppressed")
574      ENDIF
575
576      id_mbathy(:,:) = MAX( 0, id_mbathy(:,:))
577
578      ! Number of ocean level inferior or equal to jpkm1
579
580      il_ikmax = 0
581      DO jj = 1, il_jpjglo
582         DO ji = 1, il_jpiglo
583            il_ikmax = MAX( il_ikmax, id_mbathy(ji,jj) )
584         END DO
585      END DO
586
587      id_jpkmax=id_jpk
588
589      il_jpkm1=id_jpk-1
590      IF( il_ikmax > il_jpkm1 ) THEN
591         CALL logger_error("VGRID ZGR BAT CTL: maximum number of ocean level = "//&
592         &              TRIM(fct_str(il_ikmax))//" >  jpk-1."//&
593         &              " Change jpk to "//TRIM(fct_str(il_ikmax+1))//&
594         &              " to use the exact ead bathymetry" )
595      ELSE IF( il_ikmax < il_jpkm1 ) THEN
596         id_jpkmax=il_ikmax+1
597      ENDIF
598
599   END SUBROUTINE vgrid_zgr_bat_ctl
600   !-------------------------------------------------------------------
601   !> @brief This function compute bathy level in T,U,V,F point, and return
602   !> them as array of variable structure
603   !
604   !> @details
605   !> Bathymetry is read on Bathymetry file, then bathy level is computed
606   !> on T point, and finally fit to U,V,F point.
607   !>
608   !> you could specify :<br/>
609   !> - namelist where find parameter to set the depth of model levels
610   !> (default use GLORYS 75 levels parameters)
611   !> - domain structure to specify on e area to work on
612   !> - number of level to be used
613   !>
614   !> @author J.Paul
615   !> @date November, 2013 - Initial Version
616   !
617   !> @param[in] td_bathy     Bathymetry file structure
618   !> @param[in] cd_namelist  namelist
619   !> @param[in] td_dom       domain structure
620   !> @param[in] id_nlevel    number of lelvel to be used
621   !> @return array of level on T,U,V,F point (variable structure)
622   !-------------------------------------------------------------------
623   FUNCTION vgrid_get_level(td_bathy, cd_namelist, td_dom, id_nlevel)
624      IMPLICIT NONE
625      ! Argument
626      TYPE(TMPP)      , INTENT(IN) :: td_bathy
627      CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: cd_namelist
628      TYPE(TDOM)      , INTENT(IN), OPTIONAL :: td_dom
629      INTEGER(i4)     , INTENT(IN), OPTIONAL :: id_nlevel
630
631      ! function
632      TYPE(TVAR), DIMENSION(ip_npoint) :: vgrid_get_level
633
634      ! local variable
635      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_gdepw 
636      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_gdept 
637      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3w 
638      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3t
639      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3w_1d 
640      REAL(dp)   , DIMENSION(:)      , ALLOCATABLE :: dl_e3t_1d
641
642      INTEGER(i4)                                  :: il_status
643      INTEGER(i4)                                  :: il_fileid
644      INTEGER(i4)                                  :: il_jpkmax
645      INTEGER(i4), DIMENSION(2,2)                  :: il_xghost
646      INTEGER(i4), DIMENSION(:,:)    , ALLOCATABLE :: il_mbathy
647      INTEGER(i4), DIMENSION(:,:,:,:), ALLOCATABLE :: il_level
648     
649      LOGICAL                                      :: ll_exist
650
651      TYPE(TDIM) , DIMENSION(ip_maxdim)            :: tl_dim
652
653      TYPE(TDOM)                                   :: tl_dom
654
655      TYPE(TVAR)                                   :: tl_var
656
657      TYPE(TMPP)                                   :: tl_bathy
658
659      ! loop indices
660      INTEGER(i4) :: ji
661      INTEGER(i4) :: jj
662
663      INTEGER(i4) :: jip
664      INTEGER(i4) :: jjp
665
666      !namelist (intialise with GLORYS 75 levels parameters)
667      REAL(dp)                                :: dn_ppsur     = -3958.951371276829_dp
668      REAL(dp)                                :: dn_ppa0      =   103.9530096000000_dp
669      REAL(dp)                                :: dn_ppa1      =     2.4159512690000_dp
670      REAL(dp)                                :: dn_ppa2      =   100.7609285000000_dp
671      REAL(dp)                                :: dn_ppkth     =    15.3510137000000_dp
672      REAL(dp)                                :: dn_ppkth2    =    48.0298937200000_dp
673      REAL(dp)                                :: dn_ppacr     =     7.0000000000000_dp
674      REAL(dp)                                :: dn_ppacr2    =    13.000000000000_dp
675      REAL(dp)                                :: dn_ppdzmin   = 6._dp
676      REAL(dp)                                :: dn_pphmax    = 5750._dp
677      INTEGER(i4)                             :: in_nlevel    = 75
678
679      REAL(dp)                                :: dn_e3zps_min = 25._dp
680      REAL(dp)                                :: dn_e3zps_rat = 0.2_dp
681      !----------------------------------------------------------------
682      NAMELIST /namzgr/ &
683      &  dn_ppsur,     &
684      &  dn_ppa0,      &
685      &  dn_ppa1,      &
686      &  dn_ppa2,      &
687      &  dn_ppkth,     &
688      &  dn_ppkth2,    &
689      &  dn_ppacr,     &
690      &  dn_ppacr2,    &
691      &  dn_ppdzmin,   &
692      &  dn_pphmax,    &
693      &  in_nlevel
694
695      NAMELIST /namzps/ &
696      &  dn_e3zps_min, &
697      &  dn_e3zps_rat
698      !----------------------------------------------------------------
699
700      IF( PRESENT(cd_namelist) )THEN
701         !1- read namelist
702         INQUIRE(FILE=TRIM(cd_namelist), EXIST=ll_exist)
703         IF( ll_exist )THEN
704 
705            il_fileid=fct_getunit()
706
707            OPEN( il_fileid, FILE=TRIM(cd_namelist), &
708            &                FORM='FORMATTED',       &
709            &                ACCESS='SEQUENTIAL',    &
710            &                STATUS='OLD',           &
711            &                ACTION='READ',          &
712            &                IOSTAT=il_status)
713            CALL fct_err(il_status)
714            IF( il_status /= 0 )THEN
715               CALL logger_fatal("VGRID GET LEVEL: ERROR opening "//&
716               &  TRIM(cd_namelist))
717            ENDIF
718
719            READ( il_fileid, NML = namzgr )
720            READ( il_fileid, NML = namzps )
721
722            CLOSE( il_fileid, IOSTAT=il_status )
723            CALL fct_err(il_status)
724            IF( il_status /= 0 )THEN
725               CALL logger_error("VGRID GET LEVELL: ERROR closing "//&
726               &  TRIM(cd_namelist))
727            ENDIF
728
729         ELSE
730
731            CALL logger_fatal("VGRID GET LEVEL: ERROR. can not find "//&
732            &  TRIM(cd_namelist))
733
734         ENDIF
735      ENDIF
736
737      ! copy structure
738      tl_bathy=mpp_copy(td_bathy)
739
740      ! get domain
741      IF( PRESENT(td_dom) )THEN
742         tl_dom=dom_copy(td_dom)
743      ELSE
744         CALL logger_debug("VGRID GET LEVEL: get dom from "//&
745         &  TRIM(tl_bathy%c_name))
746         tl_dom=dom_init(tl_bathy)
747      ENDIF
748
749      ! get ghoste cell
750      il_xghost(:,:)=grid_get_ghost(tl_bathy)
751
752      ! open mpp files
753      CALL iom_dom_open(tl_bathy, tl_dom)
754
755      ! check namelist
756      IF( PRESENT(id_nlevel) ) in_nlevel=id_nlevel
757      IF( in_nlevel == 0 )THEN
758         CALL logger_fatal("VGRID GET LEVEL: number of level to be used "//&
759         &                 "is not specify. check namelist.")
760      ENDIF
761
762      ! read bathymetry
763      tl_var=iom_dom_read_var(tl_bathy,'bathymetry',tl_dom)
764      ! clean
765      CALL dom_clean(tl_dom)
766
767      ! remove ghost cell
768      CALL grid_del_ghost(tl_var, il_xghost(:,:))
769
770      ! force _FillValue (land) to be 0
771      WHERE( tl_var%d_value(:,:,1,1) == tl_var%d_fill )
772         tl_var%d_value(:,:,1,1)=0
773      END WHERE
774
775      ! clean
776      CALL iom_dom_close(tl_bathy)
777      CALL mpp_clean(tl_bathy)
778
779      ! compute vertical grid
780      ALLOCATE( dl_gdepw(in_nlevel), dl_gdept(in_nlevel) ) 
781      ALLOCATE(   dl_e3w(in_nlevel),   dl_e3t(in_nlevel) ) 
782      ALLOCATE(   dl_e3w_1d(in_nlevel),   dl_e3t_1d(in_nlevel) ) 
783      CALL vgrid_zgr_z( dl_gdepw(:), dl_gdept(:), dl_e3w(:), dl_e3t(:), &
784      &                 dl_e3w_1d, dl_e3t_1d, &
785      &                 dn_ppkth, dn_ppkth2, dn_ppacr, dn_ppacr2,       &
786      &                 dn_ppdzmin, dn_pphmax, &
787      &                 dn_ppa0, dn_ppa1, dn_ppa2, dn_ppsur )
788
789      ! compute bathy level on T point
790      ALLOCATE( il_mbathy(tl_var%t_dim(1)%i_len, &
791      &                   tl_var%t_dim(2)%i_len ) )
792      CALL vgrid_zgr_zps( il_mbathy(:,:), tl_var%d_value(:,:,1,1), il_jpkmax, &
793      &                   dl_gdepw(:), dl_e3t(:),               &
794      &                   dn_e3zps_min, dn_e3zps_rat )
795
796      DEALLOCATE( dl_gdepw, dl_gdept ) 
797      DEALLOCATE(   dl_e3w,   dl_e3t )
798
799      ! compute bathy level in T,U,V,F point
800      ALLOCATE( il_level(tl_var%t_dim(1)%i_len, &
801      &                  tl_var%t_dim(2)%i_len, &
802      &                  ip_npoint,1) )
803
804      DO jj=1,tl_var%t_dim(2)%i_len
805         DO ji= 1,tl_var%t_dim(1)%i_len
806
807         jip=MIN(ji+1,tl_var%t_dim(1)%i_len)
808         jjp=MIN(jj+1,tl_var%t_dim(2)%i_len)
809
810         ! T point
811         il_level(ji,jj,jp_T,1)=il_mbathy(ji,jj)
812         ! U point
813         il_level(ji,jj,jp_U,1)=MIN( il_mbathy(ji, jj ), il_mbathy(jip, jj ))
814         ! V point
815         il_level(ji,jj,jp_V,1)=MIN( il_mbathy(ji, jj ), il_mbathy(ji , jjp))
816         ! F point
817         il_level(ji,jj,jp_F,1)=MIN( il_mbathy(ji, jj ), il_mbathy(jip, jj ), &
818         &                           il_mbathy(ji, jjp), il_mbathy(jip, jjp))
819
820         ENDDO
821      ENDDO
822
823      DEALLOCATE( il_mbathy )
824
825      tl_dim(:)=dim_copy(tl_var%t_dim(:))
826      ! clean
827      CALL var_clean(tl_var)
828
829      ! only 2 first dimension to be used
830      tl_dim(3:4)%l_use=.FALSE.
831
832      vgrid_get_level(jp_T)=var_init( 'tlevel', il_level(:,:,jp_T:jp_T,:), &
833      &                                td_dim=tl_dim(:) )
834      vgrid_get_level(jp_U)=var_init( 'ulevel', il_level(:,:,jp_U:jp_U,:), &
835      &                                td_dim=tl_dim(:))
836      vgrid_get_level(jp_V)=var_init( 'vlevel', il_level(:,:,jp_V:jp_V,:), &
837      &                                td_dim=tl_dim(:))
838      vgrid_get_level(jp_F)=var_init( 'flevel', il_level(:,:,jp_F:jp_F,:), &
839      &                                td_dim=tl_dim(:))
840
841      DEALLOCATE( il_level )
842
843      CALL grid_add_ghost( vgrid_get_level(jp_T), il_xghost(:,:) )
844      CALL grid_add_ghost( vgrid_get_level(jp_U), il_xghost(:,:) )
845      CALL grid_add_ghost( vgrid_get_level(jp_V), il_xghost(:,:) )
846      CALL grid_add_ghost( vgrid_get_level(jp_F), il_xghost(:,:) )
847
848      ! clean
849      CALL dim_clean(tl_dim(:))
850
851   END FUNCTION vgrid_get_level
852END MODULE vgrid
853
Note: See TracBrowser for help on using the repository browser.