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 utils/tools/NESTING/src – NEMO

source: utils/tools/NESTING/src/agrif_readwrite.f90 @ 10093

Last change on this file since 10093 was 10093, checked in by clem, 5 years ago

just cosmetics for a better understanding of the namelist parameters
  • Property svn:keywords set to Id
File size: 39.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
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('mbathy',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('mbathy' ,dimnames,name,Grid%bathy_level,'float')
244    !
245    CALL copy_ncdf_att('nav_lon',TRIM(parent_bathy_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
246    CALL copy_ncdf_att('nav_lat',TRIM(parent_bathy_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
247    CALL copy_ncdf_att('mbathy' ,TRIM(parent_bathy_file),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_meter(name,CoarseGrid,ChildGrid)
259  !*****************************************************
260  !
261  INTEGER FUNCTION Read_Bathy_meter(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_batmet_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_meter = 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_meter = 1
414    RETURN
415    !     
416  END FUNCTION Read_Bathy_meter
417  !
418  !
419  !*****************************************************
420  !   function Read_Bathy_meter(name,CoarseGrid,ChildGrid)
421  !*****************************************************
422  !
423  INTEGER FUNCTION Read_Bathymeter(name,Grid)
424    !
425    !
426    USE io_netcdf
427    !     
428    CHARACTER(*) name
429    TYPE(Coordinates) :: Grid
430    !
431    CALL read_ncdf_var(parent_batmet_name,name,Grid%Bathy_meter)   
432    !
433    WRITE(*,*) ' '
434    WRITE(*,*) 'Reading bathymetry file: ',name
435    WRITE(*,*) ' '     
436    !
437    Read_Bathymeter = 1
438    !     
439  END FUNCTION Read_Bathymeter
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_batmet_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_batmet_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
503     REAL*8 , ALLOCATABLE, DIMENSION(:,:)   ::   ff_t, ff_f
504     INTEGER, ALLOCATABLE, DIMENSION(:,:)   ::   bottom_level, top_level
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) )
513     ALLOCATE( bottom_level(nx,ny), top_level(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     ! top/bottom levels
530     bottom_level(:,:) = Grid%bathy_level(:,:)
531     top_level(:,:)    = MIN( 1, bottom_level(:,:) )
532
533     ! vertical scale factors
534     CALL zgr_z( e3t_1d, e3w_1d, gdept_1d )   
535     DO jk = 1, N
536        e3t_0  (:,:,jk) = e3t_1d  (jk)
537        e3u_0  (:,:,jk) = e3t_1d  (jk)
538        e3v_0  (:,:,jk) = e3t_1d  (jk)
539        e3f_0  (:,:,jk) = e3t_1d  (jk)
540        e3w_0  (:,:,jk) = e3w_1d  (jk)
541        e3uw_0 (:,:,jk) = e3w_1d  (jk)
542        e3vw_0 (:,:,jk) = e3w_1d  (jk)
543     END DO
544
545     ! logicals and others
546     ln_sco = 0
547     ln_isfcav = 0
548     IF( partial_steps ) THEN
549        ln_zps = 1
550        ln_zco = 0
551     ELSE
552        ln_zps = 0
553        ln_zco = 1
554     ENDIF
555
556     ! closed domain (agrif)
557     jperio = 0
558
559     bottom_level(1:nx,1   ) = 0
560     bottom_level(1:nx,  ny) = 0
561     bottom_level(1   ,1:ny) = 0
562     bottom_level(  nx,1:ny) = 0
563
564     top_level(1:nx,1   ) = 0
565     top_level(1:nx,  ny) = 0
566     top_level(1   ,1:ny) = 0
567     top_level(  nx,1:ny) = 0
568
569     ! -------------------
570     ! write domain_cfg.nc
571     ! -------------------
572     status = nf90_create(name,NF90_WRITE,ncid)
573     status = nf90_close(ncid)
574     !
575     ! dimensions
576     dimnames = (/'x','y','z'/)
577     CALL write_ncdf_dim(dimnames(1),name,nx)
578     CALL write_ncdf_dim(dimnames(2),name,ny)
579     CALL write_ncdf_dim(dimnames(3),name,N)
580     !     
581     ! variables
582     CALL write_ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float')     
583     CALL write_ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float')
584     CALL write_ncdf_var('nav_lev',dimnames(3)  ,name,gdept_1d    ,'float')
585     !
586     CALL write_ncdf_var('jpiglo',name,nx    ,'integer')
587     CALL write_ncdf_var('jpjglo',name,ny    ,'integer')
588     CALL write_ncdf_var('jpkglo',name,N     ,'integer')
589     CALL write_ncdf_var('jperio',name,jperio,'integer')
590     !
591     CALL write_ncdf_var('ln_zco'   ,name,ln_zco   ,'integer')
592     CALL write_ncdf_var('ln_zps'   ,name,ln_zps   ,'integer')
593     CALL write_ncdf_var('ln_sco'   ,name,ln_sco   ,'integer')
594     CALL write_ncdf_var('ln_isfcav',name,ln_isfcav,'integer')
595
596     CALL write_ncdf_var('glamt',dimnames(1:2),name,Grid%glamt,'double')
597     CALL write_ncdf_var('glamu',dimnames(1:2),name,Grid%glamu,'double')
598     CALL write_ncdf_var('glamv',dimnames(1:2),name,Grid%glamv,'double')
599     CALL write_ncdf_var('glamf',dimnames(1:2),name,Grid%glamf,'double')
600     CALL write_ncdf_var('gphit',dimnames(1:2),name,Grid%gphit,'double')
601     CALL write_ncdf_var('gphiu',dimnames(1:2),name,Grid%gphiu,'double')
602     CALL write_ncdf_var('gphiv',dimnames(1:2),name,Grid%gphiv,'double')
603     CALL write_ncdf_var('gphif',dimnames(1:2),name,Grid%gphif,'double')     
604
605     CALL write_ncdf_var('e1t',dimnames(1:2),name,Grid%e1t,'double')     
606     CALL write_ncdf_var('e1u',dimnames(1:2),name,Grid%e1u,'double')     
607     CALL write_ncdf_var('e1v',dimnames(1:2),name,Grid%e1v,'double')     
608     CALL write_ncdf_var('e1f',dimnames(1:2),name,Grid%e1f,'double')
609     CALL write_ncdf_var('e2t',dimnames(1:2),name,Grid%e2t,'double')
610     CALL write_ncdf_var('e2u',dimnames(1:2),name,Grid%e2u,'double')
611     CALL write_ncdf_var('e2v',dimnames(1:2),name,Grid%e2v,'double')
612     CALL write_ncdf_var('e2f',dimnames(1:2),name,Grid%e2f,'double')
613
614     CALL write_ncdf_var('ff_f',dimnames(1:2),name,ff_f,'double')
615     CALL write_ncdf_var('ff_t',dimnames(1:2),name,ff_t,'double')
616
617     CALL write_ncdf_var('e3t_1d',dimnames(3),name,e3t_1d,'double')
618     CALL write_ncdf_var('e3w_1d',dimnames(3),name,e3w_1d,'double')
619
620     CALL write_ncdf_var('e3t_0' ,dimnames(:),name,e3t_0 ,'double')
621     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double')
622     CALL write_ncdf_var('e3u_0' ,dimnames(:),name,e3u_0 ,'double')
623     CALL write_ncdf_var('e3v_0' ,dimnames(:),name,e3v_0 ,'double')
624     CALL write_ncdf_var('e3f_0' ,dimnames(:),name,e3f_0 ,'double')
625     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double')
626     CALL write_ncdf_var('e3uw_0',dimnames(:),name,e3uw_0,'double')
627     CALL write_ncdf_var('e3vw_0',dimnames(:),name,e3vw_0,'double')
628
629     CALL write_ncdf_var('bottom_level',dimnames(1:2),name,bottom_level,'integer')
630     CALL write_ncdf_var('top_level'   ,dimnames(1:2),name,top_level   ,'integer')
631
632     CALL write_ncdf_var('bathy_meter' ,dimnames(1:2),name,Grid%bathy_meter,'float')
633
634     ! some attributes     
635     CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
636     CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
637
638     CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name)
639     CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name)
640     CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name)
641     CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name)
642     CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name)
643     CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name)
644     CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name)
645     CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name)
646
647     CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name)
648     CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name)
649     CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name)
650     CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name)
651     CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name)
652     CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name)
653     CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name)
654     CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name)           
655     !
656     ! control print
657     WRITE(*,*) ' '
658     WRITE(*,*) 'Writing domcfg file: ',name
659     WRITE(*,*) ' '
660     !
661     DEALLOCATE( ff_t, ff_f )
662     DEALLOCATE( bottom_level, top_level )
663     DEALLOCATE( e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )
664     !     
665     write_domcfg = 1
666
667  END FUNCTION write_domcfg
668  !
669  SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d)
670     !!----------------------------------------------------------------------
671     !!                   ***  ROUTINE zgr_z (from NEMO4) ***
672     !!                   
673     !! ** Purpose :   set the depth of model levels and the resulting
674     !!      vertical scale factors.
675     !!
676     !! ** Method  :   z-coordinate system (use in all type of coordinate)
677     !!        The depth of model levels is defined from an analytical
678     !!      function the derivative of which gives the scale factors.
679     !!        both depth and scale factors only depend on k (1d arrays).
680     !!              w-level: gdepw_1d  = gdep(k)
681     !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
682     !!              t-level: gdept_1d  = gdep(k+0.5)
683     !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
684     !!
685     !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
686     !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
687     !!
688     !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
689     !!----------------------------------------------------------------------
690     INTEGER  ::   jk                   ! dummy loop indices
691     INTEGER  ::   jpk
692     REAL*8 ::   zt, zw                 ! temporary scalars
693     REAL*8 ::   zsur, za0, za1, zkth   ! Values set from parameters in
694     REAL*8 ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
695     REAL*8 ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
696
697     REAL*8, DIMENSION(:), INTENT(inout) ::  e3t_1d, e3w_1d, gdept_1d
698     REAL*8, DIMENSION(N)                ::  gdepw_1d
699     !!----------------------------------------------------------------------
700     !
701     !
702     ! Set variables from parameters
703     ! ------------------------------
704     zkth = ppkth       ;   zacr = ppacr
705     zdzmin = ppdzmin   ;   zhmax = pphmax
706     zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
707
708     ! If pa1 and pa0 and psur are et to pp_to_be_computed
709     !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
710    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
711         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN
712        !
713        za1  = (  ppdzmin - pphmax / FLOAT(N-1)  )                                                      &
714           & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(N-1) * (  LOG( COSH( (N - ppkth) / ppacr) )      &
715           &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
716        za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
717        zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
718     ELSE
719        za1 = pa1 ;       za0 = pa0 ;          zsur = psur
720        za2 = pa2                            ! optional (ldbletanh=T) double tanh parameter
721     ENDIF
722
723     ! Reference z-coordinate (depth - scale factor at T- and W-points)
724     ! ======================
725     IF( ppkth == 0. ) THEN            !  uniform vertical grid
726
727        za1 = zhmax / FLOAT(N-1)
728
729        DO jk = 1, N
730           zw = FLOAT( jk )
731           zt = FLOAT( jk ) + 0.5
732           gdepw_1d(jk) = ( zw - 1 ) * za1
733           gdept_1d(jk) = ( zt - 1 ) * za1
734           e3w_1d  (jk) =  za1
735           e3t_1d  (jk) =  za1
736        END DO
737     ELSE                                ! Madec & Imbard 1996 function
738        IF( .NOT. ldbletanh ) THEN
739           DO jk = 1, N
740              zw = REAL( jk )
741              zt = REAL( jk ) + 0.5
742              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
743              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
744              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
745              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
746           END DO
747        ELSE
748           DO jk = 1, N
749              zw = FLOAT( jk )
750              zt = FLOAT( jk ) + 0.5
751              ! Double tanh function
752              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
753                 &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
754              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
755                 &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
756              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
757                 &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
758              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
759                 &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
760           END DO
761        ENDIF
762        gdepw_1d(1) = 0.                    ! force first w-level to be exactly at zero
763     ENDIF
764
765     IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
766        !
767        !==>>>   need to be like this to compute the pressure gradient with ISF.
768        !        If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
769        !        define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
770        !
771        DO jk = 1, N-1
772           e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
773        END DO
774        e3t_1d(N) = e3t_1d(N-1)   ! we don't care because this level is masked in NEMO
775
776        DO jk = 2, N
777           e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
778        END DO
779        e3w_1d(1  ) = 2. * (gdept_1d(1) - gdepw_1d(1))
780     END IF
781
782     !
783  END SUBROUTINE zgr_z
784
785
786  !*****************************************************
787  !   function set_child_name(Parentname,Childname)
788  !*****************************************************
789  !
790  SUBROUTINE set_child_name(Parentname,Childname)
791    !
792    CHARACTER(*),INTENT(in) :: Parentname
793    CHARACTER(*),INTENT(out) :: Childname
794    CHARACTER(2) :: prefix
795    INTEGER :: pos
796    !   
797    pos  = INDEX(TRIM(Parentname),'/',back=.TRUE.)
798    !
799    prefix=Parentname(pos+1:pos+2)
800    IF (prefix == '1_') THEN
801       Childname = '2_'//Parentname(pos+3:LEN(Parentname)) 
802    ELSEIF (prefix == '2_') THEN
803       Childname = '3_'//Parentname(pos+3:LEN(Parentname)) 
804    ELSEIF (prefix == '3_') THEN
805       Childname = '4_'//Parentname(pos+3:LEN(Parentname)) 
806    ELSEIF (prefix == '4_') THEN
807       Childname = '5_'//Parentname(pos+3:LEN(Parentname)) 
808    ELSE
809       Childname = '1_'//Parentname(pos+1:LEN(Parentname)) 
810    ENDIF
811    !   
812  END SUBROUTINE set_child_name
813  !
814  !*****************************************************
815  !   subroutine get_interptype(varname,interp_type,conservation)
816  !*****************************************************
817  !
818  SUBROUTINE get_interptype( varname,interp_type,conservation )
819    !
820    LOGICAL,OPTIONAL :: conservation
821    CHARACTER(*) :: interp_type,varname
822    INTEGER :: pos,pos1,pos2,pos3,i,k 
823    LOGICAL :: find       
824    i=1
825    DO WHILE ( TRIM(VAR_INTERP(i)) .NE. 'NULL' )     
826       pos = INDEX( TRIM(VAR_INTERP(i)) , TRIM(varname) )
827       IF ( pos .NE. 0 ) THEN     
828          pos1 = INDEX( TRIM(VAR_INTERP(i)) , 'bicubic' )
829          pos2 = INDEX( TRIM(VAR_INTERP(i)) , 'bilinear' )
830          pos3 = INDEX( TRIM(VAR_INTERP(i)) , 'conservative' )
831          ! initialize interp_type
832          IF( pos1 .NE. 0 ) interp_type = 'bicubic'
833          IF( pos2 .NE. 0 ) interp_type = 'bilinear'
834          IF( pos1 .EQ. 0 .AND. pos2 .EQ. 0) interp_type = 'bicubic'
835          ! initialize conservation                                           
836          IF( pos3 .NE. 0 .AND. PRESENT(conservation) ) THEN
837             conservation = .TRUE.
838             RETURN
839          ELSE
840             conservation = .FALSE.
841          ENDIF
842          find = .FALSE.
843          IF( PRESENT(conservation) ) THEN
844             k=0
845             conservation = .FALSE.         
846             DO WHILE( k < SIZE(flxtab) .AND. .NOT.find ) 
847                k = k+1
848                IF( TRIM(varname) .EQ. TRIM(flxtab(k)) ) THEN       
849                   conservation = .TRUE.
850                   find = .TRUE.
851                ENDIF
852             END DO
853          ENDIF
854          RETURN
855       ENDIF
856       i = i+1
857    END DO
858    !default values interp_type = bicubic // conservation = false     
859    interp_type = 'bicubic'     
860    IF( PRESENT(conservation) ) conservation = .FALSE.
861
862    RETURN
863    !   
864  END SUBROUTINE get_interptype
865  !     
866  !*****************************************************
867  !   subroutine Init_mask(name,Grid)
868  !*****************************************************
869  !
870  SUBROUTINE Init_mask(name,Grid,jpiglo,jpjglo)
871    !
872    USE io_netcdf
873    !     
874    CHARACTER(*) name
875    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
876    TYPE(Coordinates) :: Grid
877    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
878    !
879    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
880       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
881    ELSE
882       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
883    ENDIF
884
885
886    !
887    WRITE(*,*) 'Init masks in T,U,V,F points'   
888    !
889    nx = SIZE(Grid%Bathy_level,1)
890    ny = SIZE(Grid%Bathy_level,2)
891    !
892    !     
893    ALLOCATE(Grid%tmask(nx,ny,N), &
894         Grid%umask(nx,ny,N), &
895         Grid%vmask(nx,ny,N), &
896         Grid%fmask(nx,ny,N))
897    !
898    DO k = 1,N
899       !     
900       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
901          Grid%tmask(:,:,k) = 0
902       ELSEWHERE
903          Grid%tmask(:,:,k) = 1
904       END WHERE
905       !
906    END DO
907    !
908    Grid%umask(1:nx-1,:,:) = Grid%tmask(1:nx-1,:,:)*Grid%tmask(2:nx,:,:)
909    Grid%vmask(:,1:ny-1,:) = Grid%tmask(:,1:ny-1,:)*Grid%tmask(:,2:ny,:)
910    !
911    Grid%umask(nx,:,:) = Grid%tmask(nx,:,:)
912    Grid%vmask(:,ny,:) = Grid%tmask(:,ny,:)
913    !     
914    Grid%fmask(1:nx-1,1:ny-1,:) = Grid%tmask(1:nx-1,1:ny-1,:)*Grid%tmask(2:nx,1:ny-1,:)* &
915         Grid%tmask(1:nx-1,2:ny,:)*Grid%tmask(2:nx,2:ny,:) 
916    !
917    Grid%fmask(nx,:,:) = Grid%tmask(nx,:,:)
918    Grid%fmask(:,ny,:) = Grid%tmask(:,ny,:)
919    !
920    ALLOCATE(zwf(nx,ny))
921    !     
922    DO k = 1,N
923       !
924       zwf(:,:) = Grid%fmask(:,:,k)
925       !         
926       DO j = 2, ny-1
927          DO i = 2,nx-1           
928             IF( Grid%fmask(i,j,k) == 0. ) THEN               
929                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)))
930             END IF
931          END DO
932       END DO
933       !
934       DO j = 2, ny-1
935          IF( Grid%fmask(1,j,k) == 0. ) THEN
936             Grid%fmask(1,j,k) = shlat * MIN(1.,MAX(zwf(2,j),zwf(1,j+1),zwf(1,j-1)))
937          ENDIF
938
939          IF( Grid%fmask(nx,j,k) == 0. ) THEN
940             Grid%fmask(nx,j,k) = shlat * MIN(1.,MAX(zwf(nx,j+1),zwf(nx-1,j),zwf(nx,j-1)))
941          ENDIF
942       END DO
943       !         
944       DO i = 2, nx-1         
945          IF( Grid%fmask(i,1,k) == 0. ) THEN
946             Grid%fmask(i, 1 ,k) = shlat*MIN( 1.,MAX(zwf(i+1,1),zwf(i,2),zwf(i-1,1)))
947          ENDIF
948          !           
949          IF( Grid%fmask(i,ny,k) == 0. ) THEN
950             Grid%fmask(i,ny,k) = shlat * MIN(1.,MAX(zwf(i+1,ny),zwf(i-1,ny),zwf(i,ny-1)))
951          ENDIF
952       END DO
953       !!
954    END DO
955    !!     
956  END SUBROUTINE Init_mask
957  !
958  !*****************************************************
959  !   subroutine Init_Tmask(name,Grid)
960  !*****************************************************
961  !
962  SUBROUTINE Init_Tmask(name,Grid,jpiglo,jpjglo)
963    !
964    USE io_netcdf
965    !     
966    CHARACTER(*) name
967    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
968    TYPE(Coordinates) :: Grid
969    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
970    !
971    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
972       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
973    ELSE
974       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
975    ENDIF
976    !
977    nx = SIZE(Grid%Bathy_level,1)
978    ny = SIZE(Grid%Bathy_level,2) 
979    !
980    WRITE(*,*) 'Init masks in T points'   
981    !     
982    ALLOCATE(Grid%tmask(nx,ny,N))
983    !
984    DO k = 1,N
985       !     
986       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
987          Grid%tmask(:,:,k) = 0.
988       ELSEWHERE
989          Grid%tmask(:,:,k) = 1.
990       END WHERE
991       !
992    END DO
993    !     
994  END SUBROUTINE Init_Tmask
995  !
996  !*****************************************************
997  !   subroutine get_mask(name,Grid)
998  !*****************************************************
999  !
1000  SUBROUTINE get_mask(level,posvar,mask,filename)
1001    !
1002    USE io_netcdf
1003    !     
1004    CHARACTER(*) filename
1005    CHARACTER(*) posvar
1006    INTEGER :: level, nx, ny
1007    LOGICAL,DIMENSION(:,:),POINTER :: mask
1008    INTEGER,DIMENSION(:,:),POINTER :: maskT,maskU,maskV
1009    !     
1010    TYPE(Coordinates) :: Grid
1011    !
1012    CALL read_ncdf_var('Bathy_level',filename,Grid%Bathy_level)
1013    !
1014    nx = SIZE(Grid%Bathy_level,1)
1015    ny = SIZE(Grid%Bathy_level,2)
1016    ALLOCATE(maskT(nx,ny),mask(nx,ny))
1017    mask = .TRUE.
1018    !
1019    WHERE(Grid%Bathy_level(:,:) <= level-1 )   
1020       maskT(:,:) = 0
1021    ELSEWHERE
1022       maskT(:,:) = 1
1023    END WHERE
1024    !
1025    SELECT CASE(posvar)
1026       !
1027    CASE('T')
1028       !
1029       WHERE(maskT > 0)
1030          mask = .TRUE.
1031       ELSEWHERE
1032          mask = .FALSE.
1033       END WHERE
1034       DEALLOCATE(maskT)
1035       !
1036    CASE('U')
1037       !
1038       ALLOCATE(maskU(nx,ny))
1039       maskU(1:nx-1,:) = maskT(1:nx-1,:)*maskT(2:nx,:)
1040       maskU(nx,:) = maskT(nx,:)
1041       WHERE(maskU > 0)
1042          mask = .TRUE.
1043       ELSEWHERE
1044          mask = .FALSE.
1045       END WHERE
1046       DEALLOCATE(maskU,maskT)
1047       !
1048    CASE('V')
1049       !
1050       ALLOCATE(maskV(nx,ny))
1051       maskV(:,1:ny-1) = maskT(:,1:ny-1)*maskT(:,2:ny)
1052       maskV(:,ny) = maskT(:,ny)
1053       WHERE(maskT > 0)
1054          mask = .TRUE.
1055       ELSEWHERE
1056          mask = .FALSE.
1057       END WHERE
1058       DEALLOCATE(maskV,maskT)
1059       !
1060    END SELECT
1061    !     
1062  END SUBROUTINE get_mask
1063  !
1064  !
1065  !*****************************************************
1066  !   subroutine read_dimg_var(unit,irec,field)
1067  !*****************************************************
1068  !
1069  SUBROUTINE read_dimg_var(unit,irec,field,jpk)
1070    !     
1071    INTEGER :: unit,irec,jpk
1072    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
1073    INTEGER :: k
1074    !     
1075    DO k = 1,jpk
1076       READ(unit,REC=irec) field(:,:,k,1)
1077       irec = irec + 1
1078    ENDDO
1079    !
1080  END SUBROUTINE read_dimg_var
1081  !
1082  !
1083  !*****************************************************
1084  !   subroutine read_dimg_var(unit,irec,field)
1085  !*****************************************************
1086  !
1087  SUBROUTINE write_dimg_var(unit,irec,field,jpk)
1088    !     
1089    INTEGER :: unit,irec,jpk
1090    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
1091    INTEGER :: k
1092    !     
1093    DO k = 1,jpk
1094       WRITE(unit,REC=irec) field(:,:,k,1)
1095       irec = irec + 1
1096    ENDDO
1097    !
1098  END SUBROUTINE write_dimg_var
1099
1100END MODULE agrif_readwrite
Note: See TracBrowser for help on using the repository browser.