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 @ 10027

Last change on this file since 10027 was 10025, checked in by clem, 6 years ago

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