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

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

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