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/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src – NEMO

source: NEMO/branches/2018/dev_r5003_MERCATOR6_CRS/NEMOGCM/TOOLS/SIREN/src/vgrid.f90 @ 10115

Last change on this file since 10115 was 10115, checked in by cbricaud, 6 years ago

phase 3.6 coarsening branch with nemo_3.6_rev9192

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