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

Last change on this file since 11231 was 11231, checked in by clem, 15 months ago

remove the imposed closed boundaries when using the nesting tools

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