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/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src – NEMO

Context Navigation

source: branches/2013/dev_rev4119_MERCATOR4_CONFMAN/NEMOGCM/TOOLS/SIREN/src/vgrid.f90 @ 4213

Last change on this file since 4213 was 4213, checked in by cbricaud, 10 years ago
first draft of the CONFIGURATION MANAGER demonstrator
File size: 26.4 KB

Line
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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: