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

source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/TOOLS/SIREN/src/vgrid.f90 @ 5213

Last change on this file since 5213 was 5213, checked in by davestorkey, 9 years ago

Merge in trunk changes up to rev 5107.

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