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.
agrif_readwrite.f90 in branches/DEV_r1879_FCM/NEMOGCM/TOOLS/NESTING/src – NEMO

source: branches/DEV_r1879_FCM/NEMOGCM/TOOLS/NESTING/src/agrif_readwrite.f90 @ 2143

Last change on this file since 2143 was 2143, checked in by rblod, 12 years ago

Improvement of FCM branch

  • Property svn:keywords set to Id
File size: 32.6 KB
Line 
1!************************************************************************
2! Fortran 95 OPA Nesting tools                  *
3!                          *
4!     Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr)   *
5!                          *
6!************************************************************************
7!
8MODULE agrif_readwrite
9  !
10  USE agrif_types
11  !
12  IMPLICIT NONE
13  !
14CONTAINS
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,z
142    REAL*8,DIMENSION(:),POINTER :: tabtemp
143    INTEGER,DIMENSION(:),POINTER :: tabint
144    CHARACTER(len=20),DIMENSION(4) :: dimnames
145    !
146    status = nf90_create(name,NF90_WRITE,ncid)
147    status = nf90_close(ncid)
148    !           
149    CALL Write_Ncdf_dim('x',name,nxfin)
150    CALL Write_Ncdf_dim('y',name,nyfin)
151    IF(.NOT. iom_activated) CALL Write_Ncdf_dim('z',name,1)
152    CALL Write_Ncdf_dim('time',name,0)
153    !     
154    dimnames(1)='x'
155    dimnames(2)='y'
156    CALL Write_Ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float')     
157    CALL Write_Ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float')
158    !
159    IF(.NOT. iom_activated) THEN
160       ! copy nav_lev variable -> IOIPSL
161       CALL Read_Ncdf_dim('z',parent_coordinate_file,z)
162       ALLOCATE(tabtemp(z))
163       CALL Read_Ncdf_var('nav_lev',TRIM(parent_coordinate_file),tabtemp)
164       CALL Write_Ncdf_var('nav_lev','z',name,tabtemp,'float')           
165       DEALLOCATE(tabtemp)
166    ENDIF
167    !
168    CALL Read_Ncdf_var('time',TRIM(parent_coordinate_file),tabtemp)
169    CALL Write_Ncdf_var('time','time',name,tabtemp,'float')           
170    DEALLOCATE(tabtemp)     
171    CALL Read_Ncdf_var('time_steps',TRIM(parent_coordinate_file),tabint)
172    CALL Write_Ncdf_var('time_steps','time',name,tabint) 
173    !     
174    dimnames(1)='x'
175    dimnames(2)='y'
176    IF(iom_activated) THEN
177       dimnames(3)='time'
178    ELSE
179       dimnames(3)='z'
180       dimnames(4)='time'
181    ENDIF
182
183    CALL Write_Ncdf_var('glamt',dimnames,name,Grid%glamt,3,'double')
184    CALL Write_Ncdf_var('glamu',dimnames,name,Grid%glamu,3,'double')
185    CALL Write_Ncdf_var('glamv',dimnames,name,Grid%glamv,3,'double')
186    CALL Write_Ncdf_var('glamf',dimnames,name,Grid%glamf,3,'double')
187    CALL Write_Ncdf_var('gphit',dimnames,name,Grid%gphit,3,'double')
188    CALL Write_Ncdf_var('gphiu',dimnames,name,Grid%gphiu,3,'double')
189    CALL Write_Ncdf_var('gphiv',dimnames,name,Grid%gphiv,3,'double')
190    CALL Write_Ncdf_var('gphif',dimnames,name,Grid%gphif,3,'double')     
191    CALL Write_Ncdf_var('e1t',dimnames,name,Grid%e1t,3,'double')     
192    CALL Write_Ncdf_var('e1u',dimnames,name,Grid%e1u,3,'double')     
193    CALL Write_Ncdf_var('e1v',dimnames,name,Grid%e1v,3,'double')     
194    CALL Write_Ncdf_var('e1f',dimnames,name,Grid%e1f,3,'double')
195    CALL Write_Ncdf_var('e2t',dimnames,name,Grid%e2t,3,'double')
196    CALL Write_Ncdf_var('e2u',dimnames,name,Grid%e2u,3,'double')
197    CALL Write_Ncdf_var('e2v',dimnames,name,Grid%e2v,3,'double')
198    CALL Write_Ncdf_var('e2f',dimnames,name,Grid%e2f,3,'double')
199    !     
200    CALL Copy_Ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,&
201         MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
202    CALL Copy_Ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,&
203         MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
204    CALL Copy_Ncdf_att('nav_lev',TRIM(parent_coordinate_file),name)
205    CALL Copy_Ncdf_att('time',TRIM(parent_coordinate_file),name)
206    CALL Copy_Ncdf_att('time_steps',TRIM(parent_coordinate_file),name)
207    CALL Copy_Ncdf_att('glamt',TRIM(parent_coordinate_file),name)
208    CALL Copy_Ncdf_att('glamu',TRIM(parent_coordinate_file),name)
209    CALL Copy_Ncdf_att('glamv',TRIM(parent_coordinate_file),name)
210    CALL Copy_Ncdf_att('glamf',TRIM(parent_coordinate_file),name)
211    CALL Copy_Ncdf_att('gphit',TRIM(parent_coordinate_file),name)
212    CALL Copy_Ncdf_att('gphiu',TRIM(parent_coordinate_file),name)
213    CALL Copy_Ncdf_att('gphiv',TRIM(parent_coordinate_file),name)
214    CALL Copy_Ncdf_att('gphif',TRIM(parent_coordinate_file),name)
215    CALL Copy_Ncdf_att('e1t',TRIM(parent_coordinate_file),name)
216    CALL Copy_Ncdf_att('e1u',TRIM(parent_coordinate_file),name)
217    CALL Copy_Ncdf_att('e1v',TRIM(parent_coordinate_file),name)
218    CALL Copy_Ncdf_att('e1f',TRIM(parent_coordinate_file),name)
219    CALL Copy_Ncdf_att('e2t',TRIM(parent_coordinate_file),name)
220    CALL Copy_Ncdf_att('e2u',TRIM(parent_coordinate_file),name)
221    CALL Copy_Ncdf_att('e2v',TRIM(parent_coordinate_file),name)
222    CALL Copy_Ncdf_att('e2f',TRIM(parent_coordinate_file),name)           
223    !
224    WRITE(*,*) ' '
225    WRITE(*,*) 'Writing coordinates file: ',name
226    IF(.NOT. iom_activated) WRITE(*,*) 'IOISPL format'
227    IF(iom_activated) WRITE(*,*) 'IOM format'     
228    WRITE(*,*) ' '
229    !
230    Write_Coordinates = 1
231    !     
232  END FUNCTION Write_Coordinates
233  !
234  !
235  !
236  !*****************************************************
237  !   function Read_Bathy_level(name,Grid)
238  !*****************************************************
239  !
240  INTEGER FUNCTION Read_Bathy_level(name,Grid)
241    !
242    USE io_netcdf
243    !     
244    CHARACTER(*) name
245    TYPE(Coordinates) :: Grid
246    !
247    CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level)   
248    !
249    WRITE(*,*) ' '
250    WRITE(*,*) 'Reading bathymetry file: ',name
251    WRITE(*,*) ' '     
252    !
253    Read_Bathy_level = 1
254    !     
255  END FUNCTION Read_Bathy_level
256  !
257  !*****************************************************
258  !   function Write_Bathy_level(name,Grid)
259  !*****************************************************
260  !
261  INTEGER FUNCTION Write_Bathy_level(name,Grid)
262    !
263    USE io_netcdf
264    !     
265    CHARACTER(*) name
266    TYPE(Coordinates) :: Grid
267    INTEGER :: status,ncid,deptht
268    REAL*8,DIMENSION(:),POINTER :: timedepth_temp
269    REAL*8,DIMENSION(:,:,:),POINTER :: bathy3d
270    REAL*8,DIMENSION(:,:,:,:),POINTER :: bathy4d
271    CHARACTER(len=20),DIMENSION(4) :: dimnames
272    !
273    status = nf90_create(name,NF90_WRITE,ncid)
274    status = nf90_close(ncid)         
275    !
276    CALL Write_Ncdf_dim('x',name,nxfin)
277    CALL Write_Ncdf_dim('y',name,nyfin)
278    IF(.NOT. iom_activated ) CALL Write_Ncdf_dim('deptht',name,1)
279    CALL Write_Ncdf_dim('time_counter',name,0)
280    !     
281    dimnames(1)='x'
282    dimnames(2)='y' 
283    IF(iom_activated ) THEN
284       dimnames(3)='time_counter'
285    ELSE       
286       dimnames(3)='deptht'
287       dimnames(4)='time_counter'       
288    ENDIF
289    !     
290    CALL Write_Ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float')
291    CALL Write_Ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float')
292    !
293    ! deptht must be considered for IOIPSL format
294    !
295    IF ( Dims_Existence( 'deptht' , TRIM(parent_bathy_level) ) &
296         .AND. .NOT. iom_activated) THEN
297       CALL Read_Ncdf_dim('deptht',TRIM(parent_bathy_level),deptht)
298       ALLOCATE(timedepth_temp(deptht))
299       CALL Read_Ncdf_var('deptht',TRIM(parent_bathy_level),timedepth_temp)
300       CALL Write_Ncdf_var('deptht','deptht',name,timedepth_temp,'float')
301       CALL Copy_Ncdf_att('deptht',TRIM(parent_bathy_level),name)
302       DEALLOCATE(timedepth_temp)
303    ENDIF
304    !
305    CALL Read_Ncdf_var('time_counter',TRIM(parent_bathy_level),timedepth_temp)
306    CALL Write_Ncdf_var('time_counter','time_counter',name,timedepth_temp,'float')                 
307    !
308    IF(iom_activated) THEN
309       ALLOCATE(bathy3d(nxfin,nyfin,1))
310       bathy3d(:,:,1) = Grid%bathy_level(:,:)
311       CALL Write_Ncdf_var('Bathy_level',dimnames(1:3),name,bathy3d,'float')
312       DEALLOCATE(bathy3d)
313    ELSE
314       ALLOCATE(bathy4d(nxfin,nyfin,1,1))
315       bathy4d(:,:,1,1) = Grid%bathy_level(:,:)
316       CALL Write_Ncdf_var('Bathy_level',dimnames,name,bathy4d,'float')
317       DEALLOCATE(bathy4d)             
318    ENDIF
319    !
320    CALL Copy_Ncdf_att('nav_lon',TRIM(parent_bathy_level),name,   &
321         MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
322    CALL Copy_Ncdf_att('nav_lat',TRIM(parent_bathy_level),name,   &
323         MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
324    IF(.NOT. iom_activated )  &
325         CALL Copy_Ncdf_att('deptht',TRIM(parent_bathy_level),name)
326    CALL Copy_Ncdf_att('time_counter',TRIM(parent_bathy_level),name)
327    CALL Copy_Ncdf_att('Bathy_level',TRIM(parent_bathy_level),name)       
328    !
329    WRITE(*,*) ' '
330    WRITE(*,*) 'Writing bathymetry file: ',name
331    IF(.NOT. iom_activated) WRITE(*,*) 'IOISPL format'
332    IF(iom_activated) WRITE(*,*) 'IOM format' 
333    WRITE(*,*) ' '
334    !
335    Write_Bathy_level = 1
336    !
337  END FUNCTION Write_Bathy_level
338  !
339  !*****************************************************
340  !   function Read_Bathy_meter(name,CoarseGrid,ChildGrid)
341  !*****************************************************
342  !
343  INTEGER FUNCTION Read_Bathy_meter(name,CoarseGrid,ChildGrid,Pacifique)
344    !
345    USE io_netcdf
346    CHARACTER(*) name
347    INTEGER :: i,j,tabdim1,tabdim2
348    INTEGER, DIMENSION(1) :: i_min,i_max,j_min,j_max
349    REAL*8,POINTER,DIMENSION(:) :: topo_lon,topo_lat
350    INTEGER :: status,ncid,varid 
351    LOGICAL,OPTIONAL :: Pacifique
352    TYPE(Coordinates) :: CoarseGrid,ChildGrid
353    !
354    IF( Dims_Existence('lon',name) .AND. & 
355         Dims_Existence('lat',name)) THEN
356       WRITE(*,*) '****'
357       WRITE(*,*) ' etopo format for external high resolution database  '
358       WRITE(*,*) '****'
359       CALL Read_Ncdf_var('lon',name,topo_lon)
360       CALL Read_Ncdf_var('lat',name,topo_lat)
361    ELSE IF( Dims_Existence('x',name) &
362         .AND. Dims_Existence('y',name) ) THEN
363       WRITE(*,*) '****'
364       WRITE(*,*) ' OPA format for external high resolution database  '
365       WRITE(*,*) '****'
366       CALL Read_Ncdf_var('nav_lon',name,CoarseGrid%nav_lon)
367       CALL Read_Ncdf_var('nav_lat',name,CoarseGrid%nav_lat)
368       CALL Read_Ncdf_var('Bathymetry',name,CoarseGrid%Bathy_meter)
369       !           
370       IF ( PRESENT(Pacifique) ) THEN
371          IF(Pacifique) THEN
372             WHERE(CoarseGrid%nav_lon < 0.001) 
373                CoarseGrid%nav_lon = CoarseGrid%nav_lon + 360.
374             END WHERE
375          ENDIF
376       ENDIF
377       !     
378       Read_Bathy_meter = 1
379       RETURN     
380    ELSE
381       WRITE(*,*) '****'
382       WRITE(*,*) '*** ERROR Bad format for external high resolution database'
383       WRITE(*,*) '****'
384       STOP 
385    ENDIF
386    !
387    IF( MAXVAL(ChildGrid%glamt) > 180 ) THEN                 
388       !         
389       WHERE( topo_lon < 0 )
390          topo_lon = topo_lon + 360.
391       END WHERE
392       !         
393       i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon))
394
395       !
396       i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon))                   
397       !
398       j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat))
399       j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat))
400       !         
401       tabdim1 = ( SIZE(topo_lon) - i_min(1) + 1 ) + i_max(1)                   
402       !         
403       IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
404          tabdim2 = j_max(1)-j_min(1)+5 
405       ELSE     
406          tabdim2 = j_max(1)-j_min(1)+1   
407       ENDIF
408       !
409       ALLOCATE(CoarseGrid%nav_lon(tabdim1,tabdim2))
410       ALLOCATE(CoarseGrid%nav_lat(tabdim1,tabdim2))
411       ALLOCATE(CoarseGrid%Bathy_meter(tabdim1,tabdim2))         
412       !
413       DO j = 1,SIZE(CoarseGrid%nav_lat,1) 
414          IF( j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
415             CoarseGrid%nav_lat(j,1:tabdim2) = topo_lat(j_min(1)-2:j_max(1)+3) 
416          ELSE
417             CoarseGrid%nav_lat(j,1:tabdim2) = topo_lat(j_min(1):j_max(1)+1)
418          ENDIF
419       END DO
420       !
421       DO i = 1,SIZE(CoarseGrid%nav_lon,2) 
422          !                     
423          CoarseGrid%nav_lon(1:SIZE(topo_lon)-i_min(1),i) = topo_lon(i_min(1):SIZE(topo_lon))
424          CoarseGrid%nav_lon(1+SIZE(topo_lon)-i_min(1):tabdim1,i) = topo_lon(1:i_max(1))
425          !   
426       END DO
427       status = nf90_open(name,NF90_NOWRITE,ncid)
428       status = nf90_inq_varid(ncid,'topo',varid)
429       !         
430       status=nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter(1:SIZE(topo_lon)-i_min(1),:), &
431            start=(/i_min(1),j_min(1)-2/),count=(/SIZE(topo_lon)-i_min(1),tabdim2/))
432
433       status=nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter(1+SIZE(topo_lon)-i_min(1):tabdim1,:), &
434            start=(/1,j_min(1)-2/),count=(/i_max(1)+1,tabdim2/))               
435       !
436    ELSE
437       !
438       i_min = MAXLOC(topo_lon,mask = topo_lon < MINVAL(ChildGrid%nav_lon))
439       i_max = MINLOC(topo_lon,mask = topo_lon > MAXVAL(ChildGrid%nav_lon))
440       j_min = MAXLOC(topo_lat,mask = topo_lat < MINVAL(ChildGrid%nav_lat))
441       j_max = MINLOC(topo_lat,mask = topo_lat > MAXVAL(ChildGrid%nav_lat))
442       !     
443       IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN
444          tabdim1 = i_max(1)-i_min(1)+5 
445       ELSE     
446          tabdim1 = i_max(1)-i_min(1)+1   
447       ENDIF
448       !
449       IF(j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
450          tabdim2 = j_max(1)-j_min(1)+5 
451       ELSE     
452          tabdim2 = j_max(1)-j_min(1)+1   
453       ENDIF
454       !
455       WRITE(*,*) ' '
456       WRITE(*,*) 'Reading bathymetry file: ',name
457       WRITE(*,*) ' '
458       !
459       ALLOCATE(CoarseGrid%nav_lon(tabdim1,tabdim2))
460       ALLOCATE(CoarseGrid%nav_lat(tabdim1,tabdim2))
461       ALLOCATE(CoarseGrid%Bathy_meter(tabdim1,tabdim2))
462       !
463       DO i = 1,SIZE(CoarseGrid%nav_lon,2)
464          IF(i_min(1)-2 >= 1 .AND. i_max(1)+3 <= SIZE(topo_lon,1) ) THEN
465             CoarseGrid%nav_lon(1:tabdim1,i) = topo_lon(i_min(1)-2:i_max(1)+3)
466          ELSE     
467             CoarseGrid%nav_lon(1:tabdim1,i) = topo_lon(i_min(1):i_max(1)+1)
468          ENDIF
469          !           
470       END DO
471       !     
472       DO j = 1,SIZE(CoarseGrid%nav_lat,1) 
473          IF( j_min(1)-2 >= 1 .AND. j_max(1)+3 <= SIZE(topo_lat,1) ) THEN
474             CoarseGrid%nav_lat(j,1:tabdim2) = topo_lat(j_min(1)-2:j_max(1)+3) 
475          ELSE
476             CoarseGrid%nav_lat(j,1:tabdim2) = topo_lat(j_min(1):j_max(1)+1)
477          ENDIF
478       END DO
479       !
480       !     
481       status = nf90_open(name,NF90_NOWRITE,ncid)
482       status = nf90_inq_varid(ncid,'topo',varid)
483       IF( j_min(1)-2 >= 1 .AND. i_min(1)-2 >= 1 ) THEN
484          status=nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter, &
485               start=(/i_min(1)-2,j_min(1)-2/),count=(/tabdim1,tabdim2/))
486       ELSE
487          status=nf90_get_var(ncid,varid,CoarseGrid%Bathy_meter, &
488               start=(/i_min(1),j_min(1)/),count=(/tabdim1,tabdim2/)) 
489       ENDIF
490       !
491    ENDIF
492    !
493    !
494    status = nf90_close(ncid)     
495    !
496    WHERE(CoarseGrid%Bathy_meter.GE.0)
497       CoarseGrid%Bathy_meter = 0.0
498    END WHERE
499    !     
500    DO j = 1,SIZE(CoarseGrid%Bathy_meter,2)
501       DO i = 1,SIZE(CoarseGrid%Bathy_meter,1)
502
503          CoarseGrid%Bathy_meter(i,j) = -1.0*CoarseGrid%Bathy_meter(i,j)     
504
505       END DO
506    END DO
507    !     
508    Read_Bathy_meter = 1
509    RETURN
510    !     
511  END FUNCTION Read_Bathy_meter
512  !
513  !
514  !*****************************************************
515  !   function Read_Bathy_meter(name,CoarseGrid,ChildGrid)
516  !*****************************************************
517  !
518  INTEGER FUNCTION Read_Bathymeter(name,Grid)
519    !
520    !
521    USE io_netcdf
522    !     
523    CHARACTER(*) name
524    TYPE(Coordinates) :: Grid
525    !
526    CALL Read_Ncdf_var('Bathymetry',name,Grid%Bathy_meter)   
527    !
528    WRITE(*,*) ' '
529    WRITE(*,*) 'Reading bathymetry file: ',name
530    WRITE(*,*) ' '     
531    !
532    Read_Bathymeter = 1
533    !     
534  END FUNCTION Read_Bathymeter
535  !     
536  !*****************************************************
537  !   function Write_Bathy_meter(name,Grid)
538  !*****************************************************
539  !
540  INTEGER FUNCTION Write_Bathy_meter(name,Grid)
541    !
542    USE io_netcdf
543    !
544    CHARACTER(*) name
545    TYPE(Coordinates) :: Grid
546    INTEGER :: status,ncid,deptht
547    REAL*8,DIMENSION(:),POINTER :: timedepth_temp
548    REAL*8,DIMENSION(:,:,:),POINTER :: bathy3d
549    CHARACTER(len=20),DIMENSION(4) :: dimnames
550    REAL*8,DIMENSION(:,:,:,:),POINTER :: bathy4d
551    INTEGER :: nx,ny
552    !
553    status = nf90_create(name,NF90_WRITE,ncid)
554    status = nf90_close(ncid)     
555    !
556    nx = SIZE(Grid%bathy_meter,1)
557    ny = SIZE(Grid%bathy_meter,2)
558
559    CALL Write_Ncdf_dim('x',name,nx)
560    CALL Write_Ncdf_dim('y',name,ny)
561    IF(.NOT. iom_activated) CALL Write_Ncdf_dim('deptht',name,1)
562    IF(Dims_Existence('time_counter' , TRIM(parent_bathy_meter))) &
563         CALL Write_Ncdf_dim('time_counter',name,0)
564    !
565    dimnames(1)='x'
566    dimnames(2)='y' 
567    IF( iom_activated ) THEN
568       dimnames(3)='time_counter'
569    ELSE       
570       dimnames(3)='deptht'
571       dimnames(4)='time_counter'       
572    ENDIF
573    !     
574    CALL Write_Ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float')
575    CALL Write_Ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float')
576    !
577    IF ( Dims_Existence( 'deptht' , TRIM(parent_bathy_meter) ) &
578         .AND. .NOT. iom_activated ) THEN
579       CALL Read_Ncdf_dim('deptht',TRIM(parent_bathy_level),deptht)
580       ALLOCATE(timedepth_temp(deptht))     
581       CALL Read_Ncdf_var('deptht',TRIM(parent_bathy_meter),timedepth_temp)
582       CALL Write_Ncdf_var('deptht','deptht',name,timedepth_temp,'float')
583       CALL Copy_Ncdf_att('deptht',TRIM(parent_bathy_meter),name)
584       DEALLOCATE(timedepth_temp)
585    ENDIF
586    !
587    IF ( Dims_Existence( 'time_counter' , TRIM(parent_bathy_meter) ) ) THEN
588       CALL Read_Ncdf_var('time_counter',TRIM(parent_bathy_meter),timedepth_temp)
589       CALL Write_Ncdf_var('time_counter','time_counter',name,timedepth_temp,'float')
590       CALL Copy_Ncdf_att('time_counter',TRIM(parent_bathy_meter),name)
591       DEALLOCATE(timedepth_temp)           
592    ENDIF
593    !
594    IF(iom_activated) THEN
595       ALLOCATE(bathy3d(nx,ny,1))
596       bathy3d(:,:,1) = Grid%bathy_meter(:,:)
597       CALL Write_Ncdf_var('Bathymetry',dimnames(1:3),name,bathy3d,'float')
598       DEALLOCATE(bathy3d)
599    ELSE IF( .NOT. Dims_Existence( 'time_counter' , TRIM(parent_bathy_meter) ) .AND. &
600         .NOT. Dims_Existence( 'deptht' , TRIM(parent_bathy_meter) ) ) THEN
601       CALL Write_Ncdf_var('Bathymetry',dimnames(1:2),name,Grid%bathy_meter,'float')
602    ELSE
603       ALLOCATE(bathy4d(nx,ny,1,1))
604       bathy4d(:,:,1,1) = Grid%bathy_meter(:,:)
605       CALL Write_Ncdf_var('Bathymetry',dimnames,name,bathy4d,'float')
606       DEALLOCATE(bathy4d)     
607    ENDIF
608    !
609    CALL Copy_Ncdf_att('nav_lon',TRIM(parent_bathy_meter),name,   &
610         MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
611    CALL Copy_Ncdf_att('nav_lat',TRIM(parent_bathy_meter),name,   &
612         MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
613    !     
614    !           
615    CALL Copy_Ncdf_att('Bathymetry',TRIM(parent_bathy_meter),name)
616    !
617    WRITE(*,*) ' '
618    WRITE(*,*) 'Writing bathymetry file: ',name
619    IF(.NOT. iom_activated) WRITE(*,*) 'IOISPL format'
620    IF(iom_activated) WRITE(*,*) 'IOM format'     
621    WRITE(*,*) ' '
622    !
623    Write_Bathy_meter = 1
624    !     
625  END FUNCTION Write_Bathy_meter
626  !
627  !*****************************************************
628  !   function set_child_name(Parentname,Childname)
629  !*****************************************************
630  !
631  SUBROUTINE set_child_name(Parentname,Childname)
632    !
633    CHARACTER(*),INTENT(in) :: Parentname
634    CHARACTER(*),INTENT(out) :: Childname
635    CHARACTER(2) :: prefix
636    INTEGER :: pos
637    !   
638    pos  = INDEX(TRIM(Parentname),'/',back=.TRUE.)
639    !
640    prefix=Parentname(pos+1:pos+2)
641    IF (prefix == '1_') THEN
642       Childname = '2_'//Parentname(pos+3:LEN(Parentname)) 
643    ELSEIF (prefix == '2_') THEN
644       Childname = '3_'//Parentname(pos+3:LEN(Parentname)) 
645    ELSEIF (prefix == '3_') THEN
646       Childname = '4_'//Parentname(pos+3:LEN(Parentname)) 
647    ELSEIF (prefix == '4_') THEN
648       Childname = '5_'//Parentname(pos+3:LEN(Parentname)) 
649    ELSE
650       Childname = '1_'//Parentname(pos+1:LEN(Parentname)) 
651    ENDIF
652    !   
653  END SUBROUTINE set_child_name
654  !
655  !*****************************************************
656  !   function set_child_name(Parentname,Childname)
657  !*****************************************************
658  !
659  !*****************************************************
660  !   subroutine get_interptype(varname,interp_type,conservation)
661  !*****************************************************
662  !
663  SUBROUTINE get_interptype( varname,interp_type,conservation )
664    !
665    LOGICAL,OPTIONAL :: conservation
666    CHARACTER(*) :: interp_type,varname
667    INTEGER :: pos,pos1,pos2,pos3,i,k 
668    LOGICAL :: find       
669    i=1
670    DO WHILE ( TRIM(VAR_INTERP(i)) .NE. 'NULL' )     
671       pos = INDEX( TRIM(VAR_INTERP(i)) , TRIM(varname) )
672       IF ( pos .NE. 0 ) THEN     
673          pos1 = INDEX( TRIM(VAR_INTERP(i)) , 'bicubic' )
674          pos2 = INDEX( TRIM(VAR_INTERP(i)) , 'bilinear' )
675          pos3 = INDEX( TRIM(VAR_INTERP(i)) , 'conservative' )
676          ! initialize interp_type
677          IF( pos1 .NE. 0 ) interp_type = 'bicubic'
678          IF( pos2 .NE. 0 ) interp_type = 'bilinear'
679          IF( pos1 .EQ. 0 .AND. pos2 .EQ. 0) interp_type = 'bicubic'
680          ! initialize conservation                                           
681          IF( pos3 .NE. 0 .AND. PRESENT(conservation) ) THEN
682             conservation = .TRUE.
683             RETURN
684          ELSE
685             conservation = .FALSE.
686          ENDIF
687          find = .FALSE.
688          IF( PRESENT(conservation) ) THEN
689             k=0
690             conservation = .FALSE.         
691             DO WHILE( k <= SIZE(flxtab) .AND. .NOT.find ) 
692                k = k+1
693                IF( TRIM(varname) .EQ. TRIM(flxtab(k)) ) THEN       
694                   conservation = .TRUE.
695                   find = .TRUE.
696                ENDIF
697             END DO
698          ENDIF
699          RETURN
700       ENDIF
701       i = i+1
702    END DO
703    !default values interp_type = bicubic // conservation = false     
704    interp_type = 'bicubic'     
705    IF( PRESENT(conservation) ) conservation = .FALSE.
706
707    RETURN
708    !   
709  END SUBROUTINE get_interptype
710  !     
711  !*****************************************************             
712  !   end subroutine get_interptype
713  !*****************************************************
714  !           
715  !*****************************************************
716  !   subroutine Init_mask(name,Grid)
717  !*****************************************************
718  !
719  SUBROUTINE Init_mask(name,Grid,jpiglo,jpjglo)
720    !
721    USE io_netcdf
722    !     
723    CHARACTER(*) name
724    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
725    TYPE(Coordinates) :: Grid
726    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
727    !
728    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
729       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level)
730    ELSE
731       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
732    ENDIF
733
734
735    !
736    WRITE(*,*) 'Init masks in T,U,V,F points'   
737    !
738    nx = SIZE(Grid%Bathy_level,1)
739    ny = SIZE(Grid%Bathy_level,2)
740    !
741    !     
742    ALLOCATE(Grid%tmask(nx,ny,N), &
743         Grid%umask(nx,ny,N), &
744         Grid%vmask(nx,ny,N), &
745         Grid%fmask(nx,ny,N))
746    !
747    DO k = 1,N
748       !     
749       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
750          Grid%tmask(:,:,k) = 0
751       ELSEWHERE
752          Grid%tmask(:,:,k) = 1
753       END WHERE
754       !
755    END DO
756    !
757    Grid%umask(1:nx-1,:,:) = Grid%tmask(1:nx-1,:,:)*Grid%tmask(2:nx,:,:)
758    Grid%vmask(:,1:ny-1,:) = Grid%tmask(:,1:ny-1,:)*Grid%tmask(:,2:ny,:)
759    !
760    Grid%umask(nx,:,:) = Grid%tmask(nx,:,:)
761    Grid%vmask(:,ny,:) = Grid%tmask(:,ny,:)
762    !     
763    Grid%fmask(1:nx-1,1:ny-1,:) = Grid%tmask(1:nx-1,1:ny-1,:)*Grid%tmask(2:nx,1:ny-1,:)* &
764         Grid%tmask(1:nx-1,2:ny,:)*Grid%tmask(2:nx,2:ny,:) 
765    !
766    Grid%fmask(nx,:,:) = Grid%tmask(nx,:,:)
767    Grid%fmask(:,ny,:) = Grid%tmask(:,ny,:)
768    !
769    ALLOCATE(zwf(nx,ny))
770    !     
771    DO k = 1,N
772       !
773       zwf(:,:) = Grid%fmask(:,:,k)
774       !         
775       DO j = 2, ny-1
776          DO i = 2,nx-1           
777             IF( Grid%fmask(i,j,k) == 0. ) THEN               
778                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)))
779             END IF
780          END DO
781       END DO
782       !
783       DO j = 2, ny-1
784          IF( Grid%fmask(1,j,k) == 0. ) THEN
785             Grid%fmask(1,j,k) = shlat * MIN(1.,MAX(zwf(2,j),zwf(1,j+1),zwf(1,j-1)))
786          ENDIF
787
788          IF( Grid%fmask(nx,j,k) == 0. ) THEN
789             Grid%fmask(nx,j,k) = shlat * MIN(1.,MAX(zwf(nx,j+1),zwf(nx-1,j),zwf(nx,j-1)))
790          ENDIF
791       END DO
792       !         
793       DO i = 2, nx-1         
794          IF( Grid%fmask(i,1,k) == 0. ) THEN
795             Grid%fmask(i, 1 ,k) = shlat*MIN( 1.,MAX(zwf(i+1,1),zwf(i,2),zwf(i-1,1)))
796          ENDIF
797          !           
798          IF( Grid%fmask(i,ny,k) == 0. ) THEN
799             Grid%fmask(i,ny,k) = shlat * MIN(1.,MAX(zwf(i+1,ny),zwf(i-1,ny),zwf(i,ny-1)))
800          ENDIF
801       END DO
802       !!
803    END DO
804    !!     
805  END SUBROUTINE Init_mask
806  !
807  !*****************************************************
808  !   end subroutine Init_mask
809  !*****************************************************
810  !
811  !*****************************************************
812  !   subroutine Init_Tmask(name,Grid)
813  !*****************************************************
814  !
815  SUBROUTINE Init_Tmask(name,Grid,jpiglo,jpjglo)
816    !
817    USE io_netcdf
818    !     
819    CHARACTER(*) name
820    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
821    TYPE(Coordinates) :: Grid
822    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
823    !
824    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
825       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level)
826    ELSE
827       CALL Read_Ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
828    ENDIF
829    !
830    nx = SIZE(Grid%Bathy_level,1)
831    ny = SIZE(Grid%Bathy_level,2) 
832    !
833    WRITE(*,*) 'Init masks in T points'   
834    !     
835    ALLOCATE(Grid%tmask(nx,ny,N))
836    !
837    DO k = 1,N
838       !     
839       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
840          Grid%tmask(:,:,k) = 0.
841       ELSEWHERE
842          Grid%tmask(:,:,k) = 1.
843       END WHERE
844       !
845    END DO
846    !     
847  END SUBROUTINE Init_Tmask
848  !
849  !*****************************************************
850  !   subroutine get_mask(name,Grid)
851  !*****************************************************
852  !
853  SUBROUTINE get_mask(level,posvar,mask,filename)
854    !
855    USE io_netcdf
856    !     
857    CHARACTER(*) filename
858    CHARACTER(*) posvar
859    INTEGER :: level, nx, ny
860    LOGICAL,DIMENSION(:,:),POINTER :: mask
861    INTEGER,DIMENSION(:,:),POINTER :: maskT,maskU,maskV
862    !     
863    TYPE(Coordinates) :: Grid
864    !
865    CALL Read_Ncdf_var('Bathy_level',filename,Grid%Bathy_level)
866    !
867    nx = SIZE(Grid%Bathy_level,1)
868    ny = SIZE(Grid%Bathy_level,2)
869    ALLOCATE(maskT(nx,ny),mask(nx,ny))
870    mask = .TRUE.
871    !
872    WHERE(Grid%Bathy_level(:,:) <= level-1 )   
873       maskT(:,:) = 0
874    ELSEWHERE
875       maskT(:,:) = 1
876    END WHERE
877    !
878    SELECT CASE(posvar)
879       !
880    CASE('T')
881       !
882       WHERE(maskT > 0)
883          mask = .TRUE.
884       ELSEWHERE
885          mask = .FALSE.
886       END WHERE
887       DEALLOCATE(maskT)
888       !
889    CASE('U')
890       !
891       ALLOCATE(maskU(nx,ny))
892       maskU(1:nx-1,:) = maskT(1:nx-1,:)*maskT(2:nx,:)
893       maskU(nx,:) = maskT(nx,:)
894       WHERE(maskU > 0)
895          mask = .TRUE.
896       ELSEWHERE
897          mask = .FALSE.
898       END WHERE
899       DEALLOCATE(maskU,maskT)
900       !
901    CASE('V')
902       !
903       ALLOCATE(maskV(nx,ny))
904       maskV(:,1:ny-1) = maskT(:,1:ny-1)*maskT(:,2:ny)
905       maskV(:,ny) = maskT(:,ny)
906       WHERE(maskT > 0)
907          mask = .TRUE.
908       ELSEWHERE
909          mask = .FALSE.
910       END WHERE
911       DEALLOCATE(maskV,maskT)
912       !
913    END SELECT
914    !     
915  END SUBROUTINE get_mask
916  !
917  !*****************************************************
918  !   end subroutine get_mask
919  !*****************************************************
920  !       
921  !
922  !*****************************************************
923  !   subroutine read_dimg_var(unit,irec,field)
924  !*****************************************************
925  !
926  SUBROUTINE read_dimg_var(unit,irec,field,jpk)
927    !     
928    INTEGER :: unit,irec,jpk
929    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
930    INTEGER :: k
931    !     
932    DO k = 1,jpk
933       READ(unit,REC=irec) field(:,:,k,1)
934       irec = irec + 1
935    ENDDO
936    !
937  END SUBROUTINE read_dimg_var
938  !
939  !
940  !*****************************************************
941  !   subroutine read_dimg_var(unit,irec,field)
942  !*****************************************************
943  !
944  SUBROUTINE write_dimg_var(unit,irec,field,jpk)
945    !     
946    INTEGER :: unit,irec,jpk
947    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
948    INTEGER :: k
949    !     
950    DO k = 1,jpk
951       WRITE(unit,REC=irec) field(:,:,k,1)
952       irec = irec + 1
953    ENDDO
954    !
955  END SUBROUTINE write_dimg_var
956!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!             
957!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
958END MODULE agrif_readwrite
Note: See TracBrowser for help on using the repository browser.