1 | !************************************************************************ |
---|
2 | ! Fortran 95 OPA Nesting tools * |
---|
3 | ! * |
---|
4 | ! Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr) * |
---|
5 | ! * |
---|
6 | !************************************************************************ |
---|
7 | ! |
---|
8 | MODULE agrif_readwrite |
---|
9 | ! |
---|
10 | USE agrif_types |
---|
11 | ! |
---|
12 | IMPLICIT NONE |
---|
13 | ! |
---|
14 | CONTAINS |
---|
15 | ! |
---|
16 | !************************************************************************ |
---|
17 | ! * |
---|
18 | ! MODULE AGRIF_READWRITE * |
---|
19 | ! * |
---|
20 | ! module containing subroutine used for : * |
---|
21 | ! - Coordinates files reading/writing * |
---|
22 | ! - Bathymetry files reading/writing (meter and levels) * |
---|
23 | ! - Naming of child grid files * |
---|
24 | ! * |
---|
25 | !************************************************************************ |
---|
26 | ! |
---|
27 | !***************************************************** |
---|
28 | ! function Read_Coordinates(name,Grid) |
---|
29 | !***************************************************** |
---|
30 | |
---|
31 | INTEGER FUNCTION Read_Coordinates(name,Grid,Pacifique) |
---|
32 | ! |
---|
33 | USE io_netcdf |
---|
34 | ! |
---|
35 | ! file name to open |
---|
36 | ! |
---|
37 | CHARACTER(*) name |
---|
38 | LOGICAL,OPTIONAL :: Pacifique |
---|
39 | ! |
---|
40 | TYPE(Coordinates) :: Grid |
---|
41 | ! |
---|
42 | CALL read_ncdf_var('glamt',name,Grid%glamt) |
---|
43 | CALL read_ncdf_var('glamu',name,Grid%glamu) |
---|
44 | CALL read_ncdf_var('glamv',name,Grid%glamv) |
---|
45 | CALL read_ncdf_var('glamf',name,Grid%glamf) |
---|
46 | CALL read_ncdf_var('gphit',name,Grid%gphit) |
---|
47 | CALL read_ncdf_var('gphiu',name,Grid%gphiu) |
---|
48 | CALL read_ncdf_var('gphiv',name,Grid%gphiv) |
---|
49 | CALL read_ncdf_var('gphif',name,Grid%gphif) |
---|
50 | CALL read_ncdf_var('e1t',name,Grid%e1t) |
---|
51 | CALL read_ncdf_var('e1u',name,Grid%e1u) |
---|
52 | CALL read_ncdf_var('e1v',name,Grid%e1v) |
---|
53 | CALL read_ncdf_var('e1f',name,Grid%e1f) |
---|
54 | CALL read_ncdf_var('e2t',name,Grid%e2t) |
---|
55 | CALL read_ncdf_var('e2u',name,Grid%e2u) |
---|
56 | CALL read_ncdf_var('e2v',name,Grid%e2v) |
---|
57 | CALL read_ncdf_var('e2f',name,Grid%e2f) |
---|
58 | CALL read_ncdf_var('nav_lon',name,Grid%nav_lon) |
---|
59 | CALL read_ncdf_var('nav_lat',name,Grid%nav_lat) |
---|
60 | ! |
---|
61 | IF( PRESENT(Pacifique) )THEN |
---|
62 | IF ( Grid%glamt(1,1) > Grid%glamt(nxfin,nyfin) ) THEN |
---|
63 | Pacifique = .TRUE. |
---|
64 | WHERE ( Grid%glamt < 0 ) |
---|
65 | Grid%glamt = Grid%glamt + 360. |
---|
66 | END WHERE |
---|
67 | WHERE ( Grid%glamf < 0 ) |
---|
68 | Grid%glamf = Grid%glamf + 360. |
---|
69 | END WHERE |
---|
70 | WHERE ( Grid%glamu < 0 ) |
---|
71 | Grid%glamu = Grid%glamu + 360. |
---|
72 | END WHERE |
---|
73 | WHERE ( Grid%glamv < 0 ) |
---|
74 | Grid%glamv = Grid%glamv + 360. |
---|
75 | END WHERE |
---|
76 | WHERE ( Grid%nav_lon < 0 ) |
---|
77 | Grid%nav_lon = Grid%nav_lon + 360. |
---|
78 | END WHERE |
---|
79 | ENDIF |
---|
80 | ENDIF |
---|
81 | ! |
---|
82 | WRITE(*,*) ' ' |
---|
83 | WRITE(*,*) 'Reading coordinates file: ',name |
---|
84 | WRITE(*,*) ' ' |
---|
85 | ! |
---|
86 | Read_Coordinates = 1 |
---|
87 | ! |
---|
88 | END FUNCTION Read_Coordinates |
---|
89 | |
---|
90 | !***************************************************** |
---|
91 | ! function Read_Coordinates(name,Grid) |
---|
92 | !***************************************************** |
---|
93 | |
---|
94 | INTEGER FUNCTION Read_Local_Coordinates(name,Grid,strt,cnt) |
---|
95 | ! |
---|
96 | USE io_netcdf |
---|
97 | ! |
---|
98 | ! file name to open |
---|
99 | ! |
---|
100 | CHARACTER(*) name |
---|
101 | INTEGER, DIMENSION(2) :: strt,cnt |
---|
102 | ! |
---|
103 | TYPE(Coordinates) :: Grid |
---|
104 | ! |
---|
105 | CALL read_ncdf_var('glamt',name,Grid%glamt,strt,cnt) |
---|
106 | CALL read_ncdf_var('glamu',name,Grid%glamu,strt,cnt) |
---|
107 | CALL read_ncdf_var('glamv',name,Grid%glamv,strt,cnt) |
---|
108 | CALL read_ncdf_var('glamf',name,Grid%glamf,strt,cnt) |
---|
109 | CALL read_ncdf_var('gphit',name,Grid%gphit,strt,cnt) |
---|
110 | CALL read_ncdf_var('gphiu',name,Grid%gphiu,strt,cnt) |
---|
111 | CALL read_ncdf_var('gphiv',name,Grid%gphiv,strt,cnt) |
---|
112 | CALL read_ncdf_var('gphif',name,Grid%gphif,strt,cnt) |
---|
113 | CALL read_ncdf_var('e1t',name,Grid%e1t,strt,cnt) |
---|
114 | CALL read_ncdf_var('e1u',name,Grid%e1u,strt,cnt) |
---|
115 | CALL read_ncdf_var('e1v',name,Grid%e1v,strt,cnt) |
---|
116 | CALL read_ncdf_var('e1f',name,Grid%e1f,strt,cnt) |
---|
117 | CALL read_ncdf_var('e2t',name,Grid%e2t,strt,cnt) |
---|
118 | CALL read_ncdf_var('e2u',name,Grid%e2u,strt,cnt) |
---|
119 | CALL read_ncdf_var('e2v',name,Grid%e2v,strt,cnt) |
---|
120 | CALL read_ncdf_var('e2f',name,Grid%e2f,strt,cnt) |
---|
121 | CALL read_ncdf_var('nav_lon',name,Grid%nav_lon,strt,cnt) |
---|
122 | CALL read_ncdf_var('nav_lat',name,Grid%nav_lat,strt,cnt) |
---|
123 | ! |
---|
124 | WRITE(*,*) ' ' |
---|
125 | WRITE(*,*) 'Reading coordinates file: ',name |
---|
126 | WRITE(*,*) ' ' |
---|
127 | ! |
---|
128 | Read_Local_Coordinates = 1 |
---|
129 | ! |
---|
130 | END FUNCTION Read_Local_Coordinates |
---|
131 | |
---|
132 | !***************************************************** |
---|
133 | ! function Write_Coordinates(name,Grid) |
---|
134 | !***************************************************** |
---|
135 | |
---|
136 | INTEGER FUNCTION Write_Coordinates(name,Grid) |
---|
137 | ! |
---|
138 | USE io_netcdf |
---|
139 | CHARACTER(*) name |
---|
140 | TYPE(Coordinates) :: Grid |
---|
141 | INTEGER :: status,ncid |
---|
142 | CHARACTER(len=1),DIMENSION(2) :: dimnames |
---|
143 | ! |
---|
144 | status = nf90_create(name,NF90_WRITE,ncid) |
---|
145 | status = nf90_close(ncid) |
---|
146 | ! |
---|
147 | dimnames = (/ 'x','y' /) |
---|
148 | CALL write_ncdf_dim(dimnames(1),name,nxfin) |
---|
149 | CALL write_ncdf_dim(dimnames(2),name,nyfin) |
---|
150 | ! |
---|
151 | CALL write_ncdf_var('nav_lon',dimnames,name,Grid%nav_lon,'float') |
---|
152 | CALL write_ncdf_var('nav_lat',dimnames,name,Grid%nav_lat,'float') |
---|
153 | ! |
---|
154 | CALL write_ncdf_var('glamt',dimnames,name,Grid%glamt,'double') |
---|
155 | CALL write_ncdf_var('glamu',dimnames,name,Grid%glamu,'double') |
---|
156 | CALL write_ncdf_var('glamv',dimnames,name,Grid%glamv,'double') |
---|
157 | CALL write_ncdf_var('glamf',dimnames,name,Grid%glamf,'double') |
---|
158 | CALL write_ncdf_var('gphit',dimnames,name,Grid%gphit,'double') |
---|
159 | CALL write_ncdf_var('gphiu',dimnames,name,Grid%gphiu,'double') |
---|
160 | CALL write_ncdf_var('gphiv',dimnames,name,Grid%gphiv,'double') |
---|
161 | CALL write_ncdf_var('gphif',dimnames,name,Grid%gphif,'double') |
---|
162 | CALL write_ncdf_var('e1t',dimnames,name,Grid%e1t,'double') |
---|
163 | CALL write_ncdf_var('e1u',dimnames,name,Grid%e1u,'double') |
---|
164 | CALL write_ncdf_var('e1v',dimnames,name,Grid%e1v,'double') |
---|
165 | CALL write_ncdf_var('e1f',dimnames,name,Grid%e1f,'double') |
---|
166 | CALL write_ncdf_var('e2t',dimnames,name,Grid%e2t,'double') |
---|
167 | CALL write_ncdf_var('e2u',dimnames,name,Grid%e2u,'double') |
---|
168 | CALL write_ncdf_var('e2v',dimnames,name,Grid%e2v,'double') |
---|
169 | CALL write_ncdf_var('e2f',dimnames,name,Grid%e2f,'double') |
---|
170 | ! |
---|
171 | CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) |
---|
172 | CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) |
---|
173 | CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name) |
---|
174 | CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name) |
---|
175 | CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name) |
---|
176 | CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name) |
---|
177 | CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name) |
---|
178 | CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name) |
---|
179 | CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name) |
---|
180 | CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name) |
---|
181 | CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name) |
---|
182 | CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name) |
---|
183 | CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name) |
---|
184 | CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name) |
---|
185 | CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name) |
---|
186 | CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name) |
---|
187 | CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name) |
---|
188 | CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name) |
---|
189 | ! |
---|
190 | WRITE(*,*) ' ' |
---|
191 | WRITE(*,*) 'Writing coordinates file: ',name |
---|
192 | WRITE(*,*) ' ' |
---|
193 | ! |
---|
194 | Write_Coordinates = 1 |
---|
195 | ! |
---|
196 | END FUNCTION Write_Coordinates |
---|
197 | ! |
---|
198 | ! |
---|
199 | ! |
---|
200 | !***************************************************** |
---|
201 | ! function Read_Bathy_level(name,Grid) |
---|
202 | !***************************************************** |
---|
203 | ! |
---|
204 | INTEGER FUNCTION Read_Bathy_level(name,Grid) |
---|
205 | ! |
---|
206 | USE io_netcdf |
---|
207 | ! |
---|
208 | CHARACTER(*) name |
---|
209 | TYPE(Coordinates) :: Grid |
---|
210 | ! |
---|
211 | CALL read_ncdf_var(parent_level_name,name,Grid%Bathy_level) |
---|
212 | ! |
---|
213 | WRITE(*,*) ' ' |
---|
214 | WRITE(*,*) 'Reading bathymetry file: ',name |
---|
215 | WRITE(*,*) ' ' |
---|
216 | ! |
---|
217 | Read_Bathy_level = 1 |
---|
218 | ! |
---|
219 | END FUNCTION Read_Bathy_level |
---|
220 | ! |
---|
221 | !***************************************************** |
---|
222 | ! function Write_Bathy_level(name,Grid) |
---|
223 | !***************************************************** |
---|
224 | ! |
---|
225 | INTEGER FUNCTION Write_Bathy_level(name,Grid) |
---|
226 | ! |
---|
227 | USE io_netcdf |
---|
228 | ! |
---|
229 | CHARACTER(*) name |
---|
230 | TYPE(Coordinates) :: Grid |
---|
231 | INTEGER :: status,ncid |
---|
232 | CHARACTER(len=1),DIMENSION(2) :: dimnames |
---|
233 | ! |
---|
234 | status = nf90_create(name,NF90_WRITE,ncid) |
---|
235 | status = nf90_close(ncid) |
---|
236 | ! |
---|
237 | dimnames = (/ 'x','y' /) |
---|
238 | CALL write_ncdf_dim(dimnames(1),name,nxfin) |
---|
239 | CALL write_ncdf_dim(dimnames(2),name,nyfin) |
---|
240 | ! |
---|
241 | CALL write_ncdf_var('nav_lon' ,dimnames,name,Grid%nav_lon ,'float') |
---|
242 | CALL write_ncdf_var('nav_lat' ,dimnames,name,Grid%nav_lat ,'float') |
---|
243 | CALL write_ncdf_var(parent_level_name,dimnames,name,NINT(Grid%bathy_level),'integer') |
---|
244 | ! |
---|
245 | CALL copy_ncdf_att('nav_lon' ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) |
---|
246 | CALL copy_ncdf_att('nav_lat' ,TRIM(parent_bathy_level),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) |
---|
247 | CALL copy_ncdf_att(parent_level_name,TRIM(parent_bathy_level),name) |
---|
248 | ! |
---|
249 | WRITE(*,*) ' ' |
---|
250 | WRITE(*,*) 'Writing bathymetry (levels) in: ',name |
---|
251 | WRITE(*,*) ' ' |
---|
252 | ! |
---|
253 | Write_Bathy_level = 1 |
---|
254 | ! |
---|
255 | END FUNCTION Write_Bathy_level |
---|
256 | ! |
---|
257 | !***************************************************** |
---|
258 | ! function read_bathy_coord(name,CoarseGrid,ChildGrid) |
---|
259 | !***************************************************** |
---|
260 | ! |
---|
261 | INTEGER FUNCTION read_bathy_coord(name,CoarseGrid,ChildGrid,Pacifique) |
---|
262 | ! |
---|
263 | USE io_netcdf |
---|
264 | CHARACTER(*) name |
---|
265 | INTEGER :: i,j,tabdim1,tabdim2 |
---|
266 | INTEGER, DIMENSION(1) :: i_min,i_max,j_min,j_max |
---|
267 | REAL*8,POINTER,DIMENSION(:) :: topo_lon,topo_lat |
---|
268 | INTEGER :: status,ncid,varid |
---|
269 | LOGICAL,OPTIONAL :: Pacifique |
---|
270 | TYPE(Coordinates) :: CoarseGrid,ChildGrid |
---|
271 | REAL*8 :: zdel |
---|
272 | zdel = 0.5 ! Offset in degrees to extend extraction of bathymetry data |
---|
273 | ! |
---|
274 | IF( Dims_Existence('lon',name) .AND. Dims_Existence('lat',name) ) THEN |
---|
275 | WRITE(*,*) '****' |
---|
276 | WRITE(*,*) ' etopo format for external high resolution database ' |
---|
277 | WRITE(*,*) '****' |
---|
278 | CALL read_ncdf_var('lon',name,topo_lon) |
---|
279 | CALL read_ncdf_var('lat',name,topo_lat) |
---|
280 | ELSE IF( Dims_Existence('x',name) .AND. Dims_Existence('y',name) ) THEN |
---|
281 | WRITE(*,*) '****' |
---|
282 | WRITE(*,*) ' OPA format for external high resolution database ' |
---|
283 | WRITE(*,*) '****' |
---|
284 | CALL read_ncdf_var('nav_lon',name,CoarseGrid%nav_lon) |
---|
285 | CALL read_ncdf_var('nav_lat',name,CoarseGrid%nav_lat) |
---|
286 | CALL read_ncdf_var(parent_meter_name,name,CoarseGrid%Bathy_meter) |
---|
287 | ! |
---|
288 | IF ( PRESENT(Pacifique) ) THEN |
---|
289 | IF(Pacifique) THEN |
---|
290 | WHERE(CoarseGrid%nav_lon < 0.001) |
---|
291 | CoarseGrid%nav_lon = CoarseGrid%nav_lon + 360. |
---|
292 | END WHERE |
---|
293 | ENDIF |
---|
294 | ENDIF |
---|
295 | ! |
---|
296 | read_bathy_coord = 1 |
---|
297 | RETURN |
---|
298 | ELSE |
---|
299 | WRITE(*,*) '****' |
---|
300 | WRITE(*,*) '*** ERROR Bad format for external high resolution database' |
---|
301 | WRITE(*,*) '****' |
---|
302 | STOP |
---|
303 | ENDIF |
---|
304 | ! |
---|
305 | IF( MAXVAL(ChildGrid%glamt) > 180. ) THEN |
---|
306 | ! |
---|
307 | WHERE( topo_lon < 0. ) topo_lon = topo_lon + 360. |
---|
308 | ! |
---|
309 | i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel) |
---|
310 | i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon)+zdel) |
---|
311 | j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat)-zdel) |
---|
312 | j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)+zdel) |
---|
313 | ! |
---|
314 | tabdim1 = ( SIZE(topo_lon) - i_min(1) + 1 ) + i_max(1) |
---|
315 | ! |
---|
316 | IF( ln_agrif_domain ) THEN |
---|
317 | IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN |
---|
318 | j_min(1) = j_min(1)-2 |
---|
319 | j_max(1) = j_max(1)+3 |
---|
320 | ENDIF |
---|
321 | ENDIF |
---|
322 | tabdim2 = j_max(1) - j_min(1) + 1 |
---|
323 | ! |
---|
324 | ALLOCATE(CoarseGrid%nav_lon(tabdim1,tabdim2)) |
---|
325 | ALLOCATE(CoarseGrid%nav_lat(tabdim1,tabdim2)) |
---|
326 | ALLOCATE(CoarseGrid%Bathy_meter(tabdim1,tabdim2)) |
---|
327 | ! |
---|
328 | DO i = 1,tabdim1 |
---|
329 | CoarseGrid%nav_lat(i,:) = topo_lat(j_min(1):j_max(1)) |
---|
330 | END DO |
---|
331 | ! |
---|
332 | DO j = 1, tabdim2 |
---|
333 | ! |
---|
334 | CoarseGrid%nav_lon(1:SIZE(topo_lon)-i_min(1)+1 ,j) = topo_lon(i_min(1):SIZE(topo_lon)) |
---|
335 | CoarseGrid%nav_lon(2+SIZE(topo_lon)-i_min(1):tabdim1,j) = topo_lon(1:i_max(1)) |
---|
336 | ! |
---|
337 | END DO |
---|
338 | status = nf90_open(name,NF90_NOWRITE,ncid) |
---|
339 | status = nf90_inq_varid(ncid,elevation_name,varid) |
---|
340 | ! |
---|
341 | IF (status/=nf90_noerr) THEN |
---|
342 | WRITE(*,*)"Can't find variable: ", elevation_name |
---|
343 | STOP |
---|
344 | ENDIF |
---|
345 | ! |
---|
346 | status=nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter(1:SIZE(topo_lon)-i_min(1)+1,:), & |
---|
347 | start=(/i_min(1),j_min(1)/),count=(/SIZE(topo_lon)-i_min(1),tabdim2/)) |
---|
348 | |
---|
349 | status=nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter(2+SIZE(topo_lon)-i_min(1):tabdim1,:), & |
---|
350 | start=(/1,j_min(1)/),count=(/i_max(1),tabdim2/)) |
---|
351 | ! |
---|
352 | ELSE |
---|
353 | ! |
---|
354 | WHERE( topo_lon > 180. ) topo_lon = topo_lon - 360. |
---|
355 | ! |
---|
356 | i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon)-zdel) |
---|
357 | i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon)+zdel) |
---|
358 | j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat)-zdel) |
---|
359 | j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat)+zdel) |
---|
360 | ! |
---|
361 | IF( ln_agrif_domain ) THEN |
---|
362 | IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN |
---|
363 | i_min(1) = i_min(1)-2 |
---|
364 | i_max(1) = i_max(1)+3 |
---|
365 | ENDIF |
---|
366 | ENDIF |
---|
367 | tabdim1 = i_max(1) - i_min(1) + 1 |
---|
368 | ! |
---|
369 | IF( ln_agrif_domain ) THEN |
---|
370 | IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN |
---|
371 | j_min(1) = j_min(1)-2 |
---|
372 | j_max(1) = j_max(1)+3 |
---|
373 | ENDIF |
---|
374 | ENDIF |
---|
375 | tabdim2 = j_max(1) - j_min(1) + 1 |
---|
376 | ! |
---|
377 | WRITE(*,*) ' ' |
---|
378 | WRITE(*,*) 'Reading bathymetry file: ',name |
---|
379 | WRITE(*,*) ' ' |
---|
380 | ! |
---|
381 | ALLOCATE(CoarseGrid%nav_lon(tabdim1,tabdim2)) |
---|
382 | ALLOCATE(CoarseGrid%nav_lat(tabdim1,tabdim2)) |
---|
383 | ALLOCATE(CoarseGrid%Bathy_meter(tabdim1,tabdim2)) |
---|
384 | ! |
---|
385 | DO j = 1,tabdim2 |
---|
386 | CoarseGrid%nav_lon(:,j) = topo_lon(i_min(1):i_max(1)) |
---|
387 | END DO |
---|
388 | ! |
---|
389 | DO i = 1,tabdim1 |
---|
390 | CoarseGrid%nav_lat(i,:) = topo_lat(j_min(1):j_max(1)) |
---|
391 | END DO |
---|
392 | ! |
---|
393 | status = nf90_open(name,NF90_NOWRITE,ncid) |
---|
394 | status = nf90_inq_varid(ncid,elevation_name,varid) |
---|
395 | IF (status/=nf90_noerr) THEN |
---|
396 | WRITE(*,*)"Can't find variable: ", elevation_name |
---|
397 | STOP |
---|
398 | ENDIF |
---|
399 | |
---|
400 | status = nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter, & |
---|
401 | & start=(/i_min(1),j_min(1)/),count=(/tabdim1,tabdim2/)) |
---|
402 | ! |
---|
403 | ENDIF |
---|
404 | ! |
---|
405 | status = nf90_close(ncid) |
---|
406 | ! |
---|
407 | WHERE(CoarseGrid%Bathy_meter.GE.0) |
---|
408 | CoarseGrid%Bathy_meter = 0.0 |
---|
409 | END WHERE |
---|
410 | ! |
---|
411 | CoarseGrid%Bathy_meter(:,:) = -1.0 * CoarseGrid%Bathy_meter(:,:) |
---|
412 | ! |
---|
413 | read_bathy_coord = 1 |
---|
414 | RETURN |
---|
415 | ! |
---|
416 | END FUNCTION read_bathy_coord |
---|
417 | ! |
---|
418 | ! |
---|
419 | !***************************************************** |
---|
420 | ! function read_bathy_meter(name,CoarseGrid,ChildGrid) |
---|
421 | !***************************************************** |
---|
422 | ! |
---|
423 | INTEGER FUNCTION read_bathy_meter(name,Grid) |
---|
424 | ! |
---|
425 | ! |
---|
426 | USE io_netcdf |
---|
427 | ! |
---|
428 | CHARACTER(*) name |
---|
429 | TYPE(Coordinates) :: Grid |
---|
430 | ! |
---|
431 | CALL read_ncdf_var(parent_meter_name,name,Grid%Bathy_meter) |
---|
432 | ! |
---|
433 | WRITE(*,*) ' ' |
---|
434 | WRITE(*,*) 'Reading bathymetry file: ',name |
---|
435 | WRITE(*,*) ' ' |
---|
436 | ! |
---|
437 | read_bathy_meter = 1 |
---|
438 | ! |
---|
439 | END FUNCTION read_bathy_meter |
---|
440 | ! |
---|
441 | !***************************************************** |
---|
442 | ! function Write_Bathy_meter(name,Grid) |
---|
443 | !***************************************************** |
---|
444 | ! |
---|
445 | INTEGER FUNCTION Write_Bathy_meter(name,Grid) |
---|
446 | ! |
---|
447 | USE io_netcdf |
---|
448 | ! |
---|
449 | CHARACTER(*) name |
---|
450 | TYPE(Coordinates) :: Grid |
---|
451 | INTEGER :: status,ncid |
---|
452 | CHARACTER(len=1),DIMENSION(2) :: dimnames |
---|
453 | INTEGER :: nx,ny |
---|
454 | ! |
---|
455 | status = nf90_create(name,NF90_WRITE,ncid) |
---|
456 | status = nf90_close(ncid) |
---|
457 | ! |
---|
458 | nx = SIZE(Grid%bathy_meter,1) |
---|
459 | ny = SIZE(Grid%bathy_meter,2) |
---|
460 | dimnames = (/ 'x','y' /) |
---|
461 | |
---|
462 | CALL write_ncdf_dim(dimnames(1),name,nx) |
---|
463 | CALL write_ncdf_dim(dimnames(2),name,ny) |
---|
464 | ! |
---|
465 | CALL write_ncdf_var('nav_lon' ,dimnames,name,Grid%nav_lon ,'float') |
---|
466 | CALL write_ncdf_var('nav_lat' ,dimnames,name,Grid%nav_lat ,'float') |
---|
467 | CALL write_ncdf_var(parent_meter_name,dimnames,name,Grid%bathy_meter,'float') |
---|
468 | CALL write_ncdf_var('weight' ,dimnames,name,Grid%wgt ,'float') |
---|
469 | ! |
---|
470 | CALL copy_ncdf_att('nav_lon' ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) |
---|
471 | CALL copy_ncdf_att('nav_lat' ,TRIM(parent_bathy_meter),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) |
---|
472 | CALL copy_ncdf_att(parent_meter_name,TRIM(parent_bathy_meter),name) |
---|
473 | ! |
---|
474 | WRITE(*,*) ' ' |
---|
475 | WRITE(*,*) 'Writing bathymetry (meters) in: ',name |
---|
476 | WRITE(*,*) ' ' |
---|
477 | ! |
---|
478 | Write_Bathy_meter = 1 |
---|
479 | ! |
---|
480 | END FUNCTION Write_Bathy_meter |
---|
481 | ! |
---|
482 | !***************************************************** |
---|
483 | ! function write_domcfg(name,Grid) |
---|
484 | !***************************************************** |
---|
485 | |
---|
486 | INTEGER FUNCTION write_domcfg(name,Grid) |
---|
487 | !----------------------------------------- |
---|
488 | ! It creates a domain_cfg.nc used in NEMO4 |
---|
489 | !----------------------------------------- |
---|
490 | ! |
---|
491 | USE io_netcdf |
---|
492 | ! |
---|
493 | CHARACTER(*) name |
---|
494 | TYPE(Coordinates) :: Grid |
---|
495 | ! |
---|
496 | INTEGER :: status, ncid |
---|
497 | INTEGER :: nx, ny, jk |
---|
498 | INTEGER :: ln_sco, ln_isfcav, ln_zco, ln_zps, jperio |
---|
499 | REAL*8 :: rpi, rad, rday, rsiyea, rsiday, omega |
---|
500 | ! |
---|
501 | CHARACTER(len=1), DIMENSION(3) :: dimnames |
---|
502 | REAL*8 , DIMENSION(N) :: e3t_1d, e3w_1d, gdept_1d, gdepw_1d |
---|
503 | REAL*8 , ALLOCATABLE, DIMENSION(:,:) :: ff_t, ff_f, zbathy |
---|
504 | INTEGER, ALLOCATABLE, DIMENSION(:,:) :: bottom_level, top_level, mbathy |
---|
505 | REAL*8 , ALLOCATABLE, DIMENSION(:,:,:) :: e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 |
---|
506 | ! |
---|
507 | ! size of the Grid |
---|
508 | nx = SIZE(Grid%bathy_meter,1) |
---|
509 | ny = SIZE(Grid%bathy_meter,2) |
---|
510 | |
---|
511 | ! allocate needed arrays for domain_cfg |
---|
512 | ALLOCATE( ff_t(nx,ny), ff_f(nx,ny), zbathy(nx,ny) ) |
---|
513 | ALLOCATE( bottom_level(nx,ny), top_level(nx,ny), mbathy(nx,ny) ) |
---|
514 | ALLOCATE( e3t_0(nx,ny,N), e3u_0 (nx,ny,N), e3v_0 (nx,ny,N), e3f_0(nx,ny,N), & |
---|
515 | & e3w_0(nx,ny,N), e3uw_0(nx,ny,N), e3vw_0(nx,ny,N) ) |
---|
516 | |
---|
517 | ! some physical parameters |
---|
518 | rpi = 3.141592653589793 |
---|
519 | rad = 3.141592653589793 / 180. |
---|
520 | rday = 24.*60.*60. |
---|
521 | rsiyea = 365.25 * rday * 2. * rpi / 6.283076 |
---|
522 | rsiday = rday / ( 1. + rday / rsiyea ) |
---|
523 | omega = 2. * rpi / rsiday |
---|
524 | |
---|
525 | ! Coriolis |
---|
526 | ff_f(:,:) = 2. * omega * SIN( rad * Grid%gphif(:,:) ) ! compute it on the sphere at f-point |
---|
527 | ff_t(:,:) = 2. * omega * SIN( rad * Grid%gphit(:,:) ) ! - - - at t-point |
---|
528 | |
---|
529 | ! bathy in meters |
---|
530 | zbathy(:,:) = Grid%bathy_meter |
---|
531 | |
---|
532 | ! vertical scale factors |
---|
533 | CALL zgr_z( e3t_1d, e3w_1d, gdept_1d, gdepw_1d ) |
---|
534 | ! DO jk = 1, N |
---|
535 | ! e3t_0 (:,:,jk) = e3t_1d (jk) |
---|
536 | ! e3u_0 (:,:,jk) = e3t_1d (jk) |
---|
537 | ! e3v_0 (:,:,jk) = e3t_1d (jk) |
---|
538 | ! e3f_0 (:,:,jk) = e3t_1d (jk) |
---|
539 | ! e3w_0 (:,:,jk) = e3w_1d (jk) |
---|
540 | ! e3uw_0 (:,:,jk) = e3w_1d (jk) |
---|
541 | ! e3vw_0 (:,:,jk) = e3w_1d (jk) |
---|
542 | ! END DO |
---|
543 | |
---|
544 | ! logicals and others |
---|
545 | ln_sco = 0 |
---|
546 | ln_isfcav = 0 |
---|
547 | !!$ IF( partial_steps ) THEN |
---|
548 | ln_zps = 1 |
---|
549 | ln_zco = 0 |
---|
550 | !!$ ELSE |
---|
551 | !!$ ln_zps = 0 |
---|
552 | !!$ ln_zco = 1 |
---|
553 | !!$ ENDIF |
---|
554 | |
---|
555 | CALL zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & ! input |
---|
556 | & mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) ! output |
---|
557 | |
---|
558 | ! top/bottom levels |
---|
559 | bottom_level(:,:) = mbathy(:,:) |
---|
560 | top_level(:,:) = MIN( 1, mbathy(:,:) ) |
---|
561 | |
---|
562 | ! closed domain for child grid and namelist defined domain for parent grid |
---|
563 | IF( TRIM(name) == TRIM(parent_domcfg_updated) .OR. TRIM(name) == TRIM(parent_domcfg_out) ) THEN |
---|
564 | jperio = parent_jperio |
---|
565 | ELSE |
---|
566 | jperio = 0 |
---|
567 | ENDIF |
---|
568 | |
---|
569 | !IF( .NOT.ln_agrif_domain ) THEN |
---|
570 | ! bottom_level(1:nx,1 ) = 0 |
---|
571 | ! bottom_level(1:nx, ny) = 0 |
---|
572 | ! bottom_level(1 ,1:ny) = 0 |
---|
573 | ! bottom_level( nx,1:ny) = 0 |
---|
574 | |
---|
575 | ! top_level(1:nx,1 ) = 0 |
---|
576 | ! top_level(1:nx, ny) = 0 |
---|
577 | ! top_level(1 ,1:ny) = 0 |
---|
578 | ! top_level( nx,1:ny) = 0 |
---|
579 | |
---|
580 | ! zbathy(1:nx,1 ) = 0. |
---|
581 | ! zbathy(1:nx, ny) = 0. |
---|
582 | ! zbathy(1 ,1:ny) = 0. |
---|
583 | ! zbathy( nx,1:ny) = 0. |
---|
584 | !ENDIF |
---|
585 | |
---|
586 | ! ------------------- |
---|
587 | ! write domain_cfg.nc |
---|
588 | ! ------------------- |
---|
589 | status = nf90_create(name,IOR(NF90_64BIT_OFFSET,NF90_WRITE),ncid) |
---|
590 | status = nf90_close(ncid) |
---|
591 | ! |
---|
592 | ! dimensions |
---|
593 | dimnames = (/'x','y','z'/) |
---|
594 | CALL write_ncdf_dim(dimnames(1),name,nx) |
---|
595 | CALL write_ncdf_dim(dimnames(2),name,ny) |
---|
596 | CALL write_ncdf_dim(dimnames(3),name,N) |
---|
597 | ! |
---|
598 | ! variables |
---|
599 | CALL write_ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float') |
---|
600 | CALL write_ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float') |
---|
601 | CALL write_ncdf_var('nav_lev',dimnames(3) ,name,gdept_1d ,'float') |
---|
602 | ! |
---|
603 | CALL write_ncdf_var('jpiglo',name,nx ,'integer') |
---|
604 | CALL write_ncdf_var('jpjglo',name,ny ,'integer') |
---|
605 | CALL write_ncdf_var('jpkglo',name,N ,'integer') |
---|
606 | CALL write_ncdf_var('jperio',name,jperio,'integer') |
---|
607 | ! |
---|
608 | CALL write_ncdf_var('ln_zco' ,name,ln_zco ,'integer') |
---|
609 | CALL write_ncdf_var('ln_zps' ,name,ln_zps ,'integer') |
---|
610 | CALL write_ncdf_var('ln_sco' ,name,ln_sco ,'integer') |
---|
611 | CALL write_ncdf_var('ln_isfcav',name,ln_isfcav,'integer') |
---|
612 | |
---|
613 | CALL write_ncdf_var('glamt',dimnames(1:2),name,Grid%glamt,'double') |
---|
614 | CALL write_ncdf_var('glamu',dimnames(1:2),name,Grid%glamu,'double') |
---|
615 | CALL write_ncdf_var('glamv',dimnames(1:2),name,Grid%glamv,'double') |
---|
616 | CALL write_ncdf_var('glamf',dimnames(1:2),name,Grid%glamf,'double') |
---|
617 | CALL write_ncdf_var('gphit',dimnames(1:2),name,Grid%gphit,'double') |
---|
618 | CALL write_ncdf_var('gphiu',dimnames(1:2),name,Grid%gphiu,'double') |
---|
619 | CALL write_ncdf_var('gphiv',dimnames(1:2),name,Grid%gphiv,'double') |
---|
620 | CALL write_ncdf_var('gphif',dimnames(1:2),name,Grid%gphif,'double') |
---|
621 | |
---|
622 | CALL write_ncdf_var('e1t',dimnames(1:2),name,Grid%e1t,'double') |
---|
623 | CALL write_ncdf_var('e1u',dimnames(1:2),name,Grid%e1u,'double') |
---|
624 | CALL write_ncdf_var('e1v',dimnames(1:2),name,Grid%e1v,'double') |
---|
625 | CALL write_ncdf_var('e1f',dimnames(1:2),name,Grid%e1f,'double') |
---|
626 | CALL write_ncdf_var('e2t',dimnames(1:2),name,Grid%e2t,'double') |
---|
627 | CALL write_ncdf_var('e2u',dimnames(1:2),name,Grid%e2u,'double') |
---|
628 | CALL write_ncdf_var('e2v',dimnames(1:2),name,Grid%e2v,'double') |
---|
629 | CALL write_ncdf_var('e2f',dimnames(1:2),name,Grid%e2f,'double') |
---|
630 | |
---|
631 | CALL write_ncdf_var('ff_f',dimnames(1:2),name,ff_f,'double') |
---|
632 | CALL write_ncdf_var('ff_t',dimnames(1:2),name,ff_t,'double') |
---|
633 | |
---|
634 | CALL write_ncdf_var('e3t_1d',dimnames(3),name,e3t_1d,'double') |
---|
635 | CALL write_ncdf_var('e3w_1d',dimnames(3),name,e3w_1d,'double') |
---|
636 | |
---|
637 | CALL write_ncdf_var('e3t_0' ,dimnames(:),name,e3t_0 ,'double') |
---|
638 | CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double') |
---|
639 | CALL write_ncdf_var('e3u_0' ,dimnames(:),name,e3u_0 ,'double') |
---|
640 | CALL write_ncdf_var('e3v_0' ,dimnames(:),name,e3v_0 ,'double') |
---|
641 | CALL write_ncdf_var('e3f_0' ,dimnames(:),name,e3f_0 ,'double') |
---|
642 | CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double') |
---|
643 | CALL write_ncdf_var('e3uw_0',dimnames(:),name,e3uw_0,'double') |
---|
644 | CALL write_ncdf_var('e3vw_0',dimnames(:),name,e3vw_0,'double') |
---|
645 | |
---|
646 | CALL write_ncdf_var('bottom_level',dimnames(1:2),name,bottom_level,'integer') |
---|
647 | CALL write_ncdf_var('top_level' ,dimnames(1:2),name,top_level ,'integer') |
---|
648 | CALL write_ncdf_var('bathy_meter' ,dimnames(1:2),name,zbathy ,'float') |
---|
649 | |
---|
650 | ! some attributes |
---|
651 | CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon)) |
---|
652 | CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat)) |
---|
653 | |
---|
654 | CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name) |
---|
655 | CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name) |
---|
656 | CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name) |
---|
657 | CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name) |
---|
658 | CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name) |
---|
659 | CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name) |
---|
660 | CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name) |
---|
661 | CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name) |
---|
662 | |
---|
663 | CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name) |
---|
664 | CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name) |
---|
665 | CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name) |
---|
666 | CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name) |
---|
667 | CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name) |
---|
668 | CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name) |
---|
669 | CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name) |
---|
670 | CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name) |
---|
671 | ! |
---|
672 | ! control print |
---|
673 | WRITE(*,*) ' ' |
---|
674 | WRITE(*,*) 'Writing domcfg file: ',name |
---|
675 | WRITE(*,*) ' ' |
---|
676 | ! |
---|
677 | DEALLOCATE( ff_t, ff_f, zbathy ) |
---|
678 | DEALLOCATE( bottom_level, top_level, mbathy ) |
---|
679 | DEALLOCATE( e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) |
---|
680 | ! |
---|
681 | write_domcfg = 1 |
---|
682 | |
---|
683 | END FUNCTION write_domcfg |
---|
684 | ! |
---|
685 | SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d, gdepw_1d ) |
---|
686 | !!---------------------------------------------------------------------- |
---|
687 | !! *** ROUTINE zgr_z (from NEMO4) *** |
---|
688 | !! |
---|
689 | !! ** Purpose : set the depth of model levels and the resulting |
---|
690 | !! vertical scale factors. |
---|
691 | !! |
---|
692 | !! ** Method : z-coordinate system (use in all type of coordinate) |
---|
693 | !! The depth of model levels is defined from an analytical |
---|
694 | !! function the derivative of which gives the scale factors. |
---|
695 | !! both depth and scale factors only depend on k (1d arrays). |
---|
696 | !! w-level: gdepw_1d = gdep(k) |
---|
697 | !! e3w_1d(k) = dk(gdep)(k) = e3(k) |
---|
698 | !! t-level: gdept_1d = gdep(k+0.5) |
---|
699 | !! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5) |
---|
700 | !! |
---|
701 | !! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m) |
---|
702 | !! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m) |
---|
703 | !! |
---|
704 | !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. |
---|
705 | !!---------------------------------------------------------------------- |
---|
706 | INTEGER :: jk ! dummy loop indices |
---|
707 | REAL*8 :: zt, zw ! temporary scalars |
---|
708 | REAL*8 :: zsur, za0, za1, zkth ! Values set from parameters in |
---|
709 | REAL*8 :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90 |
---|
710 | REAL*8 :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters |
---|
711 | |
---|
712 | REAL*8, DIMENSION(:), INTENT(out) :: e3t_1d, e3w_1d, gdept_1d, gdepw_1d |
---|
713 | !!---------------------------------------------------------------------- |
---|
714 | ! |
---|
715 | ! |
---|
716 | ! Set variables from parameters |
---|
717 | ! ------------------------------ |
---|
718 | zkth = ppkth ; zacr = ppacr |
---|
719 | zdzmin = ppdzmin ; zhmax = pphmax |
---|
720 | zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters |
---|
721 | |
---|
722 | ! If pa1 and pa0 and psur are et to pp_to_be_computed |
---|
723 | ! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr |
---|
724 | IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) & |
---|
725 | .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN |
---|
726 | ! |
---|
727 | za1 = ( ppdzmin - pphmax / FLOAT(N-1) ) & |
---|
728 | & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(N-1) * ( LOG( COSH( (N - ppkth) / ppacr) ) & |
---|
729 | & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) |
---|
730 | za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) |
---|
731 | zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) |
---|
732 | ELSE |
---|
733 | za1 = pa1 ; za0 = pa0 ; zsur = psur |
---|
734 | za2 = pa2 ! optional (ldbletanh=T) double tanh parameter |
---|
735 | ENDIF |
---|
736 | |
---|
737 | ! Reference z-coordinate (depth - scale factor at T- and W-points) |
---|
738 | ! ====================== |
---|
739 | IF( ppkth == 0. ) THEN ! uniform vertical grid |
---|
740 | |
---|
741 | za1 = zhmax / FLOAT(N-1) |
---|
742 | |
---|
743 | DO jk = 1, N |
---|
744 | zw = FLOAT( jk ) |
---|
745 | zt = FLOAT( jk ) + 0.5 |
---|
746 | gdepw_1d(jk) = ( zw - 1 ) * za1 |
---|
747 | gdept_1d(jk) = ( zt - 1 ) * za1 |
---|
748 | e3w_1d (jk) = za1 |
---|
749 | e3t_1d (jk) = za1 |
---|
750 | END DO |
---|
751 | ELSE ! Madec & Imbard 1996 function |
---|
752 | IF( .NOT. ldbletanh ) THEN |
---|
753 | DO jk = 1, N |
---|
754 | zw = REAL( jk ) |
---|
755 | zt = REAL( jk ) + 0.5 |
---|
756 | gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) ) |
---|
757 | gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) ) |
---|
758 | e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr ) |
---|
759 | e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr ) |
---|
760 | END DO |
---|
761 | ELSE |
---|
762 | DO jk = 1, N |
---|
763 | zw = FLOAT( jk ) |
---|
764 | zt = FLOAT( jk ) + 0.5 |
---|
765 | ! Double tanh function |
---|
766 | gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) & |
---|
767 | & + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) ) |
---|
768 | gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) & |
---|
769 | & + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) ) |
---|
770 | e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) & |
---|
771 | & + za2 * TANH( (zw-zkth2) / zacr2 ) |
---|
772 | e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) & |
---|
773 | & + za2 * TANH( (zt-zkth2) / zacr2 ) |
---|
774 | END DO |
---|
775 | ENDIF |
---|
776 | gdepw_1d(1) = 0. ! force first w-level to be exactly at zero |
---|
777 | ENDIF |
---|
778 | |
---|
779 | IF ( ln_e3_dep ) THEN ! e3. = dk[gdep] |
---|
780 | ! |
---|
781 | !==>>> need to be like this to compute the pressure gradient with ISF. |
---|
782 | ! If not, level beneath the ISF are not aligned (sum(e3t) /= depth) |
---|
783 | ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively |
---|
784 | ! |
---|
785 | DO jk = 1, N-1 |
---|
786 | e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) |
---|
787 | END DO |
---|
788 | e3t_1d(N) = e3t_1d(N-1) ! we don't care because this level is masked in NEMO |
---|
789 | |
---|
790 | DO jk = 2, N |
---|
791 | e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) |
---|
792 | END DO |
---|
793 | e3w_1d(1 ) = 2. * (gdept_1d(1) - gdepw_1d(1)) |
---|
794 | END IF |
---|
795 | |
---|
796 | ! |
---|
797 | END SUBROUTINE zgr_z |
---|
798 | |
---|
799 | SUBROUTINE zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, & |
---|
800 | & mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 ) |
---|
801 | !!---------------------------------------------------------------------- |
---|
802 | !! *** ROUTINE zgr_zps (from NEMO4) *** |
---|
803 | !! |
---|
804 | !! ** Purpose : the depth and vertical scale factor in partial step |
---|
805 | !! reference z-coordinate case |
---|
806 | !! |
---|
807 | !! ** Method : Partial steps : computes the 3D vertical scale factors |
---|
808 | !! of T-, U-, V-, W-, UW-, VW and F-points that are associated with |
---|
809 | !! a partial step representation of bottom topography. |
---|
810 | !! |
---|
811 | !! The reference depth of model levels is defined from an analytical |
---|
812 | !! function the derivative of which gives the reference vertical |
---|
813 | !! scale factors. |
---|
814 | !! From depth and scale factors reference, we compute there new value |
---|
815 | !! with partial steps on 3d arrays ( i, j, k ). |
---|
816 | !! |
---|
817 | !! w-level: gdepw_0(i,j,k) = gdep(k) |
---|
818 | !! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k) |
---|
819 | !! t-level: gdept_0(i,j,k) = gdep(k+0.5) |
---|
820 | !! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5) |
---|
821 | !! |
---|
822 | !! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), |
---|
823 | !! we find the mbathy index of the depth at each grid point. |
---|
824 | !! This leads us to three cases: |
---|
825 | !! |
---|
826 | !! - bathy = 0 => mbathy = 0 |
---|
827 | !! - 1 < mbathy < jpkm1 |
---|
828 | !! - bathy > gdepw_0(jpk) => mbathy = jpkm1 |
---|
829 | !! |
---|
830 | !! Then, for each case, we find the new depth at t- and w- levels |
---|
831 | !! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw- |
---|
832 | !! and f-points. |
---|
833 | !! |
---|
834 | !! This routine is given as an example, it must be modified |
---|
835 | !! following the user s desiderata. nevertheless, the output as |
---|
836 | !! well as the way to compute the model levels and scale factors |
---|
837 | !! must be respected in order to insure second order accuracy |
---|
838 | !! schemes. |
---|
839 | !! |
---|
840 | !! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives |
---|
841 | !! - - - - - - - gdept_0, gdepw_0 and e3. are positives |
---|
842 | !! |
---|
843 | !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. |
---|
844 | !!---------------------------------------------------------------------- |
---|
845 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
846 | INTEGER :: ik |
---|
847 | INTEGER :: ikb, ikt ! temporary integers |
---|
848 | REAL*8 :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points |
---|
849 | REAL*8 :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t |
---|
850 | REAL*8 :: zmax ! temporary scalar |
---|
851 | ! |
---|
852 | REAL*8, ALLOCATABLE, DIMENSION(:,:,:) :: gdept_0, gdepw_0, zprt |
---|
853 | ! |
---|
854 | INTEGER, INTENT(in ) :: nx, ny |
---|
855 | REAL*8 , DIMENSION(:) , INTENT(in ) :: gdept_1d, gdepw_1d, e3t_1d, e3w_1d |
---|
856 | REAL*8 , DIMENSION(:,:) , INTENT(inout) :: zbathy |
---|
857 | INTEGER, DIMENSION(:,:) , INTENT( out) :: mbathy |
---|
858 | REAL*8 , DIMENSION(:,:,:), INTENT( out) :: e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 |
---|
859 | ! |
---|
860 | !!--------------------------------------------------------------------- |
---|
861 | |
---|
862 | ALLOCATE( zprt(nx,ny,N), gdept_0(nx,ny,N), gdepw_0(nx,ny,N) ) |
---|
863 | ! |
---|
864 | ! bathymetry in level (from bathy_meter) |
---|
865 | ! =================== |
---|
866 | zmax = gdepw_1d(N) + e3t_1d(N) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) ) |
---|
867 | zbathy(:,:) = MIN( zmax , zbathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) |
---|
868 | WHERE( zbathy(:,:) == 0. ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 |
---|
869 | ELSEWHERE ; mbathy(:,:) = N-1 ! ocean : initialize mbathy to the max ocean level |
---|
870 | END WHERE |
---|
871 | |
---|
872 | !!$ IF( partial_steps ) THEN |
---|
873 | ! Compute mbathy for ocean points (i.e. the number of ocean levels) |
---|
874 | ! find the number of ocean levels such that the last level thickness |
---|
875 | ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where |
---|
876 | ! e3t_1d is the reference level thickness |
---|
877 | DO jk = N-1, 1, -1 |
---|
878 | zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) |
---|
879 | WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 |
---|
880 | END DO |
---|
881 | !!$ ELSE |
---|
882 | !!$ DO jk = 1, N |
---|
883 | !!$ WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) >= gdepw_1d(jk) ) mbathy(:,:) = jk-1 |
---|
884 | !!$ END DO |
---|
885 | !!$ ENDIF |
---|
886 | |
---|
887 | ! Scale factors and depth at T- and W-points |
---|
888 | DO jk = 1, N ! intitialization to the reference z-coordinate |
---|
889 | gdept_0(:,:,jk) = gdept_1d(jk) |
---|
890 | gdepw_0(:,:,jk) = gdepw_1d(jk) |
---|
891 | e3t_0 (:,:,jk) = e3t_1d (jk) |
---|
892 | e3w_0 (:,:,jk) = e3w_1d (jk) |
---|
893 | END DO |
---|
894 | |
---|
895 | ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf |
---|
896 | !!clem: not implemented yet |
---|
897 | !! IF ( ln_isfcav == 1 ) CALL zgr_isf |
---|
898 | ! |
---|
899 | ! Scale factors and depth at T- and W-points |
---|
900 | ! IF ( ln_isfcav == 0 ) THEN |
---|
901 | DO jj = 1, ny |
---|
902 | DO ji = 1, nx |
---|
903 | ik = mbathy(ji,jj) |
---|
904 | IF( ik > 0 ) THEN ! ocean point only |
---|
905 | ! max ocean level case |
---|
906 | IF( ik == N-1 ) THEN |
---|
907 | zdepwp = zbathy(ji,jj) |
---|
908 | ze3tp = zbathy(ji,jj) - gdepw_1d(ik) |
---|
909 | ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) ) |
---|
910 | e3t_0(ji,jj,ik ) = ze3tp |
---|
911 | e3t_0(ji,jj,ik+1) = ze3tp |
---|
912 | e3w_0(ji,jj,ik ) = ze3wp |
---|
913 | e3w_0(ji,jj,ik+1) = ze3tp |
---|
914 | gdepw_0(ji,jj,ik+1) = zdepwp |
---|
915 | gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp |
---|
916 | gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp |
---|
917 | ! |
---|
918 | ELSE ! standard case |
---|
919 | IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = zbathy(ji,jj) |
---|
920 | ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) |
---|
921 | ENDIF |
---|
922 | !gm Bug? check the gdepw_1d |
---|
923 | ! ... on ik |
---|
924 | gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & |
---|
925 | & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & |
---|
926 | & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) |
---|
927 | e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & |
---|
928 | & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) |
---|
929 | e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) ) & |
---|
930 | & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) |
---|
931 | ! ... on ik+1 |
---|
932 | e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) |
---|
933 | e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) |
---|
934 | gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) |
---|
935 | ENDIF |
---|
936 | ENDIF |
---|
937 | END DO |
---|
938 | END DO |
---|
939 | ! |
---|
940 | ! DO jj = 1, ny |
---|
941 | ! DO ji = 1, nx |
---|
942 | ! ik = mbathy(ji,jj) |
---|
943 | ! IF( ik > 0 ) THEN ! ocean point only |
---|
944 | ! e3tp (ji,jj) = e3t_0(ji,jj,ik) |
---|
945 | ! e3wp (ji,jj) = e3w_0(ji,jj,ik) |
---|
946 | ! ENDIF |
---|
947 | ! END DO |
---|
948 | ! END DO |
---|
949 | ! END IF |
---|
950 | ! |
---|
951 | ! Scale factors and depth at U-, V-, UW and VW-points |
---|
952 | DO jk = 1, N ! initialisation to z-scale factors |
---|
953 | e3u_0 (:,:,jk) = e3t_1d(jk) |
---|
954 | e3v_0 (:,:,jk) = e3t_1d(jk) |
---|
955 | e3uw_0(:,:,jk) = e3w_1d(jk) |
---|
956 | e3vw_0(:,:,jk) = e3w_1d(jk) |
---|
957 | END DO |
---|
958 | |
---|
959 | DO jk = 1,N ! Computed as the minimum of neighbooring scale factors |
---|
960 | DO jj = 1, ny-1 |
---|
961 | DO ji = 1, nx-1 ! vector opt. |
---|
962 | e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) |
---|
963 | e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) |
---|
964 | e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) |
---|
965 | e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) |
---|
966 | END DO |
---|
967 | END DO |
---|
968 | END DO |
---|
969 | |
---|
970 | ! IF ( ln_isfcav == 1 ) THEN |
---|
971 | ! ! (ISF) define e3uw (adapted for 2 cells in the water column) |
---|
972 | ! DO jj = 2, ny-1 |
---|
973 | ! DO ji = 2, nx-1 ! vector opt. |
---|
974 | ! ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) |
---|
975 | ! ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) |
---|
976 | ! IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & |
---|
977 | ! & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) |
---|
978 | ! ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) |
---|
979 | ! ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) |
---|
980 | ! IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & |
---|
981 | ! & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) |
---|
982 | ! END DO |
---|
983 | ! END DO |
---|
984 | ! END IF |
---|
985 | ! |
---|
986 | |
---|
987 | DO jk = 1, N ! set to z-scale factor if zero (i.e. along closed boundaries) |
---|
988 | WHERE( e3u_0 (:,:,jk) == 0. ) e3u_0 (:,:,jk) = e3t_1d(jk) |
---|
989 | WHERE( e3v_0 (:,:,jk) == 0. ) e3v_0 (:,:,jk) = e3t_1d(jk) |
---|
990 | WHERE( e3uw_0(:,:,jk) == 0. ) e3uw_0(:,:,jk) = e3w_1d(jk) |
---|
991 | WHERE( e3vw_0(:,:,jk) == 0. ) e3vw_0(:,:,jk) = e3w_1d(jk) |
---|
992 | END DO |
---|
993 | |
---|
994 | ! Scale factor at F-point |
---|
995 | DO jk = 1, N ! initialisation to z-scale factors |
---|
996 | e3f_0(:,:,jk) = e3t_1d(jk) |
---|
997 | END DO |
---|
998 | DO jk = 1, N ! Computed as the minimum of neighbooring V-scale factors |
---|
999 | DO jj = 1, ny-1 |
---|
1000 | DO ji = 1, nx-1 ! vector opt. |
---|
1001 | e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) |
---|
1002 | END DO |
---|
1003 | END DO |
---|
1004 | END DO |
---|
1005 | ! |
---|
1006 | DO jk = 1, N ! set to z-scale factor if zero (i.e. along closed boundaries) |
---|
1007 | WHERE( e3f_0(:,:,jk) == 0. ) e3f_0(:,:,jk) = e3t_1d(jk) |
---|
1008 | END DO |
---|
1009 | |
---|
1010 | ! Control of the sign |
---|
1011 | IF( MINVAL( e3t_0 (:,:,:) ) <= 0. ) STOP ' zgr_zps : e r r o r e3t_0 <= 0' |
---|
1012 | IF( MINVAL( e3w_0 (:,:,:) ) <= 0. ) STOP ' zgr_zps : e r r o r e3w_0 <= 0' |
---|
1013 | IF( MINVAL( gdept_0(:,:,:) ) < 0. ) STOP ' zgr_zps : e r r o r gdept_0 < 0' |
---|
1014 | IF( MINVAL( gdepw_0(:,:,:) ) < 0. ) STOP ' zgr_zps : e r r o r gdepw_0 < 0' |
---|
1015 | |
---|
1016 | ! Compute gde3w_0 (vertical sum of e3w) |
---|
1017 | ! IF ( ln_isfcav ==1 ) THEN ! if cavity |
---|
1018 | ! WHERE( misfdep == 0 ) misfdep = 1 |
---|
1019 | ! DO jj = 1,ny |
---|
1020 | ! DO ji = 1,nx |
---|
1021 | ! gde3w_0(ji,jj,1) = 0.5 * e3w_0(ji,jj,1) |
---|
1022 | ! DO jk = 2, misfdep(ji,jj) |
---|
1023 | ! gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) |
---|
1024 | ! END DO |
---|
1025 | ! IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5 * e3w_0(ji,jj,misfdep(ji,jj)) |
---|
1026 | ! DO jk = misfdep(ji,jj) + 1, N |
---|
1027 | ! gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) |
---|
1028 | ! END DO |
---|
1029 | ! END DO |
---|
1030 | ! END DO |
---|
1031 | ! ELSE ! no cavity |
---|
1032 | ! gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) |
---|
1033 | ! DO jk = 2, N |
---|
1034 | ! gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) |
---|
1035 | ! END DO |
---|
1036 | ! END IF |
---|
1037 | ! |
---|
1038 | DEALLOCATE( zprt, gdept_0, gdepw_0 ) |
---|
1039 | ! |
---|
1040 | END SUBROUTINE zgr_zps |
---|
1041 | |
---|
1042 | !***************************************************** |
---|
1043 | ! function set_child_name(Parentname,Childname) |
---|
1044 | !***************************************************** |
---|
1045 | ! |
---|
1046 | SUBROUTINE set_child_name(Parentname,Childname) |
---|
1047 | ! |
---|
1048 | CHARACTER(*),INTENT(in) :: Parentname |
---|
1049 | CHARACTER(*),INTENT(out) :: Childname |
---|
1050 | CHARACTER(2) :: prefix |
---|
1051 | INTEGER :: pos |
---|
1052 | ! |
---|
1053 | pos = INDEX(TRIM(Parentname),'/',back=.TRUE.) |
---|
1054 | ! |
---|
1055 | prefix=Parentname(pos+1:pos+2) |
---|
1056 | IF (prefix == '1_') THEN |
---|
1057 | Childname = '2_'//Parentname(pos+3:LEN(Parentname)) |
---|
1058 | ELSEIF (prefix == '2_') THEN |
---|
1059 | Childname = '3_'//Parentname(pos+3:LEN(Parentname)) |
---|
1060 | ELSEIF (prefix == '3_') THEN |
---|
1061 | Childname = '4_'//Parentname(pos+3:LEN(Parentname)) |
---|
1062 | ELSEIF (prefix == '4_') THEN |
---|
1063 | Childname = '5_'//Parentname(pos+3:LEN(Parentname)) |
---|
1064 | ELSE |
---|
1065 | Childname = '1_'//Parentname(pos+1:LEN(Parentname)) |
---|
1066 | ENDIF |
---|
1067 | ! |
---|
1068 | END SUBROUTINE set_child_name |
---|
1069 | ! |
---|
1070 | !***************************************************** |
---|
1071 | ! subroutine get_interptype(varname,interp_type,conservation) |
---|
1072 | !***************************************************** |
---|
1073 | ! |
---|
1074 | SUBROUTINE get_interptype( varname,interp_type,conservation ) |
---|
1075 | ! |
---|
1076 | LOGICAL,OPTIONAL :: conservation |
---|
1077 | CHARACTER(*) :: interp_type,varname |
---|
1078 | INTEGER :: pos,pos1,pos2,pos3,i,k |
---|
1079 | LOGICAL :: find |
---|
1080 | i=1 |
---|
1081 | DO WHILE ( TRIM(VAR_INTERP(i)) /= '' ) |
---|
1082 | pos = INDEX( TRIM(VAR_INTERP(i)) , TRIM(varname) ) |
---|
1083 | IF ( pos .NE. 0 ) THEN |
---|
1084 | pos1 = INDEX( TRIM(VAR_INTERP(i)) , 'bicubic' ) |
---|
1085 | pos2 = INDEX( TRIM(VAR_INTERP(i)) , 'bilinear' ) |
---|
1086 | pos3 = INDEX( TRIM(VAR_INTERP(i)) , 'conservative' ) |
---|
1087 | ! initialize interp_type |
---|
1088 | IF( pos1 .NE. 0 ) interp_type = 'bicubic' |
---|
1089 | IF( pos2 .NE. 0 ) interp_type = 'bilinear' |
---|
1090 | IF( pos1 .EQ. 0 .AND. pos2 .EQ. 0) interp_type = 'bicubic' |
---|
1091 | ! initialize conservation |
---|
1092 | IF( pos3 .NE. 0 .AND. PRESENT(conservation) ) THEN |
---|
1093 | conservation = .TRUE. |
---|
1094 | RETURN |
---|
1095 | ELSE |
---|
1096 | conservation = .FALSE. |
---|
1097 | ENDIF |
---|
1098 | find = .FALSE. |
---|
1099 | IF( PRESENT(conservation) ) THEN |
---|
1100 | k=0 |
---|
1101 | conservation = .FALSE. |
---|
1102 | DO WHILE( k < SIZE(flxtab) .AND. .NOT.find ) |
---|
1103 | k = k+1 |
---|
1104 | IF( TRIM(varname) .EQ. TRIM(flxtab(k)) ) THEN |
---|
1105 | conservation = .TRUE. |
---|
1106 | find = .TRUE. |
---|
1107 | ENDIF |
---|
1108 | END DO |
---|
1109 | ENDIF |
---|
1110 | RETURN |
---|
1111 | ENDIF |
---|
1112 | i = i+1 |
---|
1113 | END DO |
---|
1114 | !default values interp_type = bicubic // conservation = false |
---|
1115 | interp_type = 'bicubic' |
---|
1116 | IF( PRESENT(conservation) ) conservation = .FALSE. |
---|
1117 | |
---|
1118 | RETURN |
---|
1119 | ! |
---|
1120 | END SUBROUTINE get_interptype |
---|
1121 | ! |
---|
1122 | !***************************************************** |
---|
1123 | ! subroutine Init_mask(name,Grid) |
---|
1124 | !***************************************************** |
---|
1125 | ! |
---|
1126 | SUBROUTINE Init_mask(name,Grid,jpiglo,jpjglo) |
---|
1127 | ! |
---|
1128 | USE io_netcdf |
---|
1129 | ! |
---|
1130 | CHARACTER(*) name |
---|
1131 | INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo |
---|
1132 | TYPE(Coordinates) :: Grid |
---|
1133 | REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL() |
---|
1134 | ! |
---|
1135 | IF(jpiglo == 1 .AND. jpjglo == 1) THEN |
---|
1136 | CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level) |
---|
1137 | ELSE |
---|
1138 | CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) ) |
---|
1139 | ENDIF |
---|
1140 | |
---|
1141 | |
---|
1142 | ! |
---|
1143 | WRITE(*,*) 'Init masks in T,U,V,F points' |
---|
1144 | ! |
---|
1145 | nx = SIZE(Grid%Bathy_level,1) |
---|
1146 | ny = SIZE(Grid%Bathy_level,2) |
---|
1147 | ! |
---|
1148 | ! |
---|
1149 | ALLOCATE(Grid%tmask(nx,ny,N), & |
---|
1150 | Grid%umask(nx,ny,N), & |
---|
1151 | Grid%vmask(nx,ny,N), & |
---|
1152 | Grid%fmask(nx,ny,N)) |
---|
1153 | ! |
---|
1154 | DO k = 1,N |
---|
1155 | ! |
---|
1156 | WHERE(Grid%Bathy_level(:,:) <= k-1 ) |
---|
1157 | Grid%tmask(:,:,k) = 0 |
---|
1158 | ELSEWHERE |
---|
1159 | Grid%tmask(:,:,k) = 1 |
---|
1160 | END WHERE |
---|
1161 | ! |
---|
1162 | END DO |
---|
1163 | ! |
---|
1164 | Grid%umask(1:nx-1,:,:) = Grid%tmask(1:nx-1,:,:)*Grid%tmask(2:nx,:,:) |
---|
1165 | Grid%vmask(:,1:ny-1,:) = Grid%tmask(:,1:ny-1,:)*Grid%tmask(:,2:ny,:) |
---|
1166 | ! |
---|
1167 | Grid%umask(nx,:,:) = Grid%tmask(nx,:,:) |
---|
1168 | Grid%vmask(:,ny,:) = Grid%tmask(:,ny,:) |
---|
1169 | ! |
---|
1170 | Grid%fmask(1:nx-1,1:ny-1,:) = Grid%tmask(1:nx-1,1:ny-1,:)*Grid%tmask(2:nx,1:ny-1,:)* & |
---|
1171 | Grid%tmask(1:nx-1,2:ny,:)*Grid%tmask(2:nx,2:ny,:) |
---|
1172 | ! |
---|
1173 | Grid%fmask(nx,:,:) = Grid%tmask(nx,:,:) |
---|
1174 | Grid%fmask(:,ny,:) = Grid%tmask(:,ny,:) |
---|
1175 | ! |
---|
1176 | ALLOCATE(zwf(nx,ny)) |
---|
1177 | ! |
---|
1178 | DO k = 1,N |
---|
1179 | ! |
---|
1180 | zwf(:,:) = Grid%fmask(:,:,k) |
---|
1181 | ! |
---|
1182 | DO j = 2, ny-1 |
---|
1183 | DO i = 2,nx-1 |
---|
1184 | IF( Grid%fmask(i,j,k) == 0. ) THEN |
---|
1185 | Grid%fmask(i,j,k) = shlat * MIN(1.,MAX( zwf(i+1,j),zwf(i,j+1),zwf(i-1,j),zwf(i,j-1))) |
---|
1186 | END IF |
---|
1187 | END DO |
---|
1188 | END DO |
---|
1189 | ! |
---|
1190 | DO j = 2, ny-1 |
---|
1191 | IF( Grid%fmask(1,j,k) == 0. ) THEN |
---|
1192 | Grid%fmask(1,j,k) = shlat * MIN(1.,MAX(zwf(2,j),zwf(1,j+1),zwf(1,j-1))) |
---|
1193 | ENDIF |
---|
1194 | |
---|
1195 | IF( Grid%fmask(nx,j,k) == 0. ) THEN |
---|
1196 | Grid%fmask(nx,j,k) = shlat * MIN(1.,MAX(zwf(nx,j+1),zwf(nx-1,j),zwf(nx,j-1))) |
---|
1197 | ENDIF |
---|
1198 | END DO |
---|
1199 | ! |
---|
1200 | DO i = 2, nx-1 |
---|
1201 | IF( Grid%fmask(i,1,k) == 0. ) THEN |
---|
1202 | Grid%fmask(i, 1 ,k) = shlat*MIN( 1.,MAX(zwf(i+1,1),zwf(i,2),zwf(i-1,1))) |
---|
1203 | ENDIF |
---|
1204 | ! |
---|
1205 | IF( Grid%fmask(i,ny,k) == 0. ) THEN |
---|
1206 | Grid%fmask(i,ny,k) = shlat * MIN(1.,MAX(zwf(i+1,ny),zwf(i-1,ny),zwf(i,ny-1))) |
---|
1207 | ENDIF |
---|
1208 | END DO |
---|
1209 | !! |
---|
1210 | END DO |
---|
1211 | !! |
---|
1212 | END SUBROUTINE Init_mask |
---|
1213 | ! |
---|
1214 | !***************************************************** |
---|
1215 | ! subroutine Init_Tmask(name,Grid) |
---|
1216 | !***************************************************** |
---|
1217 | ! |
---|
1218 | SUBROUTINE Init_Tmask(name,Grid,jpiglo,jpjglo) |
---|
1219 | ! |
---|
1220 | USE io_netcdf |
---|
1221 | ! |
---|
1222 | CHARACTER(*) name |
---|
1223 | INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo |
---|
1224 | TYPE(Coordinates) :: Grid |
---|
1225 | REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL() |
---|
1226 | ! |
---|
1227 | IF(jpiglo == 1 .AND. jpjglo == 1) THEN |
---|
1228 | CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level) |
---|
1229 | ELSE |
---|
1230 | CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) ) |
---|
1231 | ENDIF |
---|
1232 | ! |
---|
1233 | nx = SIZE(Grid%Bathy_level,1) |
---|
1234 | ny = SIZE(Grid%Bathy_level,2) |
---|
1235 | ! |
---|
1236 | WRITE(*,*) 'Init masks in T points' |
---|
1237 | ! |
---|
1238 | ALLOCATE( Grid%tmask(nx,ny,N) ) |
---|
1239 | ! |
---|
1240 | DO k = 1,N |
---|
1241 | ! |
---|
1242 | WHERE(Grid%Bathy_level(:,:) <= k-1 ) |
---|
1243 | Grid%tmask(:,:,k) = 0. |
---|
1244 | ELSEWHERE |
---|
1245 | Grid%tmask(:,:,k) = 1. |
---|
1246 | END WHERE |
---|
1247 | ! |
---|
1248 | END DO |
---|
1249 | ! |
---|
1250 | END SUBROUTINE Init_Tmask |
---|
1251 | ! |
---|
1252 | !***************************************************** |
---|
1253 | ! subroutine Init_ssmask(name,Grid) |
---|
1254 | !***************************************************** |
---|
1255 | ! |
---|
1256 | ! SUBROUTINE Init_ssmask(varname,filename,Grid,jpiglo,jpjglo) |
---|
1257 | ! ! |
---|
1258 | ! USE io_netcdf |
---|
1259 | ! ! |
---|
1260 | ! CHARACTER(*) varname,filename |
---|
1261 | ! INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo |
---|
1262 | ! TYPE(Coordinates) :: Grid |
---|
1263 | ! REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL() |
---|
1264 | ! ! |
---|
1265 | ! IF(jpiglo == 1 .AND. jpjglo == 1) THEN |
---|
1266 | ! CALL read_ncdf_var(varname,filename,Grid%bathy_level) |
---|
1267 | ! ELSE |
---|
1268 | ! CALL read_ncdf_var(varname,filename,Grid%bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) ) |
---|
1269 | ! ENDIF |
---|
1270 | ! ! |
---|
1271 | ! nx = SIZE(Grid%bathy_level,1) |
---|
1272 | ! ny = SIZE(Grid%bathy_level,2) |
---|
1273 | ! ! |
---|
1274 | ! WRITE(*,*) 'Init surface masks in T points' |
---|
1275 | ! ! |
---|
1276 | ! ALLOCATE( Grid%ssmask(nx,ny), Grid%ssumask(nx,ny), Grid%ssvmask(nx,ny) ) |
---|
1277 | ! ! |
---|
1278 | ! WHERE(Grid%bathy_level(:,:) <= 0. ) |
---|
1279 | ! Grid%ssmask(:,:) = 0. |
---|
1280 | ! ELSEWHERE |
---|
1281 | ! Grid%ssmask(:,:) = 1. |
---|
1282 | ! END WHERE |
---|
1283 | ! ! |
---|
1284 | ! Grid%ssumask(1:nx-1,:) = Grid%ssmask(1:nx-1,:)*Grid%ssmask(2:nx,:) |
---|
1285 | ! Grid%ssvmask(:,1:ny-1) = Grid%ssmask(:,1:ny-1)*Grid%ssmask(:,2:ny) |
---|
1286 | ! ! |
---|
1287 | ! Grid%ssumask(nx,:) = Grid%ssmask(nx,:) |
---|
1288 | ! Grid%ssvmask(:,ny) = Grid%ssmask(:,ny) |
---|
1289 | ! ! |
---|
1290 | ! END SUBROUTINE Init_ssmask |
---|
1291 | ! |
---|
1292 | !***************************************************** |
---|
1293 | ! subroutine get_mask(name,Grid) |
---|
1294 | !***************************************************** |
---|
1295 | ! |
---|
1296 | SUBROUTINE get_mask(level,posvar,mask,filename) |
---|
1297 | ! |
---|
1298 | USE io_netcdf |
---|
1299 | ! |
---|
1300 | CHARACTER(*) filename |
---|
1301 | CHARACTER(*) posvar |
---|
1302 | INTEGER :: level, nx, ny |
---|
1303 | LOGICAL,DIMENSION(:,:),POINTER :: mask |
---|
1304 | INTEGER,DIMENSION(:,:),POINTER :: maskT,maskU,maskV |
---|
1305 | ! |
---|
1306 | TYPE(Coordinates) :: Grid |
---|
1307 | ! |
---|
1308 | CALL read_ncdf_var('Bathy_level',filename,Grid%Bathy_level) |
---|
1309 | ! |
---|
1310 | nx = SIZE(Grid%Bathy_level,1) |
---|
1311 | ny = SIZE(Grid%Bathy_level,2) |
---|
1312 | ALLOCATE(maskT(nx,ny),mask(nx,ny)) |
---|
1313 | mask = .TRUE. |
---|
1314 | ! |
---|
1315 | WHERE(Grid%Bathy_level(:,:) <= level-1 ) |
---|
1316 | maskT(:,:) = 0 |
---|
1317 | ELSEWHERE |
---|
1318 | maskT(:,:) = 1 |
---|
1319 | END WHERE |
---|
1320 | ! |
---|
1321 | SELECT CASE(posvar) |
---|
1322 | ! |
---|
1323 | CASE('T') |
---|
1324 | ! |
---|
1325 | WHERE(maskT > 0) |
---|
1326 | mask = .TRUE. |
---|
1327 | ELSEWHERE |
---|
1328 | mask = .FALSE. |
---|
1329 | END WHERE |
---|
1330 | DEALLOCATE(maskT) |
---|
1331 | ! |
---|
1332 | CASE('U') |
---|
1333 | ! |
---|
1334 | ALLOCATE(maskU(nx,ny)) |
---|
1335 | maskU(1:nx-1,:) = maskT(1:nx-1,:)*maskT(2:nx,:) |
---|
1336 | maskU(nx,:) = maskT(nx,:) |
---|
1337 | WHERE(maskU > 0) |
---|
1338 | mask = .TRUE. |
---|
1339 | ELSEWHERE |
---|
1340 | mask = .FALSE. |
---|
1341 | END WHERE |
---|
1342 | DEALLOCATE(maskU,maskT) |
---|
1343 | ! |
---|
1344 | CASE('V') |
---|
1345 | ! |
---|
1346 | ALLOCATE(maskV(nx,ny)) |
---|
1347 | maskV(:,1:ny-1) = maskT(:,1:ny-1)*maskT(:,2:ny) |
---|
1348 | maskV(:,ny) = maskT(:,ny) |
---|
1349 | WHERE(maskT > 0) |
---|
1350 | mask = .TRUE. |
---|
1351 | ELSEWHERE |
---|
1352 | mask = .FALSE. |
---|
1353 | END WHERE |
---|
1354 | DEALLOCATE(maskV,maskT) |
---|
1355 | ! |
---|
1356 | END SELECT |
---|
1357 | ! |
---|
1358 | END SUBROUTINE get_mask |
---|
1359 | ! |
---|
1360 | ! |
---|
1361 | !***************************************************** |
---|
1362 | ! subroutine read_dimg_var(unit,irec,field) |
---|
1363 | !***************************************************** |
---|
1364 | ! |
---|
1365 | SUBROUTINE read_dimg_var(unit,irec,field,jpk) |
---|
1366 | ! |
---|
1367 | INTEGER :: unit,irec,jpk |
---|
1368 | REAL*8,DIMENSION(:,:,:,:),POINTER :: field |
---|
1369 | INTEGER :: k |
---|
1370 | ! |
---|
1371 | DO k = 1,jpk |
---|
1372 | READ(unit,REC=irec) field(:,:,k,1) |
---|
1373 | irec = irec + 1 |
---|
1374 | ENDDO |
---|
1375 | ! |
---|
1376 | END SUBROUTINE read_dimg_var |
---|
1377 | ! |
---|
1378 | ! |
---|
1379 | !***************************************************** |
---|
1380 | ! subroutine read_dimg_var(unit,irec,field) |
---|
1381 | !***************************************************** |
---|
1382 | ! |
---|
1383 | SUBROUTINE write_dimg_var(unit,irec,field,jpk) |
---|
1384 | ! |
---|
1385 | INTEGER :: unit,irec,jpk |
---|
1386 | REAL*8,DIMENSION(:,:,:,:),POINTER :: field |
---|
1387 | INTEGER :: k |
---|
1388 | ! |
---|
1389 | DO k = 1,jpk |
---|
1390 | WRITE(unit,REC=irec) field(:,:,k,1) |
---|
1391 | irec = irec + 1 |
---|
1392 | ENDDO |
---|
1393 | ! |
---|
1394 | END SUBROUTINE write_dimg_var |
---|
1395 | |
---|
1396 | END MODULE agrif_readwrite |
---|