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 utils/tools/SIREN/src – NEMO

source: utils/tools/SIREN/src/vgrid.f90 @ 12080

Last change on this file since 12080 was 12080, checked in by jpaul, 4 years ago

update nemo trunk

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