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