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

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

solve ticket #2280 by adding a new namelist parameter in the NESTING namelist to account for jperio of the parent grid

  • Property svn:keywords set to Id
File size: 53.0 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     bottom_level(1:nx,1   ) = 0
569     bottom_level(1:nx,  ny) = 0
570     bottom_level(1   ,1:ny) = 0
571     bottom_level(  nx,1:ny) = 0
572
573     top_level(1:nx,1   ) = 0
574     top_level(1:nx,  ny) = 0
575     top_level(1   ,1:ny) = 0
576     top_level(  nx,1:ny) = 0
577
578     zbathy(1:nx,1   ) = 0.
579     zbathy(1:nx,  ny) = 0.
580     zbathy(1   ,1:ny) = 0.
581     zbathy(  nx,1:ny) = 0.     
582     
583     ! -------------------
584     ! write domain_cfg.nc
585     ! -------------------
586     status = nf90_create(name,IOR(NF90_64BIT_OFFSET,NF90_WRITE),ncid)
587     status = nf90_close(ncid)
588     !
589     ! dimensions
590     dimnames = (/'x','y','z'/)
591     CALL write_ncdf_dim(dimnames(1),name,nx)
592     CALL write_ncdf_dim(dimnames(2),name,ny)
593     CALL write_ncdf_dim(dimnames(3),name,N)
594     !     
595     ! variables
596     CALL write_ncdf_var('nav_lon',dimnames(1:2),name,Grid%nav_lon,'float')     
597     CALL write_ncdf_var('nav_lat',dimnames(1:2),name,Grid%nav_lat,'float')
598     CALL write_ncdf_var('nav_lev',dimnames(3)  ,name,gdept_1d    ,'float')
599     !
600     CALL write_ncdf_var('jpiglo',name,nx    ,'integer')
601     CALL write_ncdf_var('jpjglo',name,ny    ,'integer')
602     CALL write_ncdf_var('jpkglo',name,N     ,'integer')
603     CALL write_ncdf_var('jperio',name,jperio,'integer')
604     !
605     CALL write_ncdf_var('ln_zco'   ,name,ln_zco   ,'integer')
606     CALL write_ncdf_var('ln_zps'   ,name,ln_zps   ,'integer')
607     CALL write_ncdf_var('ln_sco'   ,name,ln_sco   ,'integer')
608     CALL write_ncdf_var('ln_isfcav',name,ln_isfcav,'integer')
609
610     CALL write_ncdf_var('glamt',dimnames(1:2),name,Grid%glamt,'double')
611     CALL write_ncdf_var('glamu',dimnames(1:2),name,Grid%glamu,'double')
612     CALL write_ncdf_var('glamv',dimnames(1:2),name,Grid%glamv,'double')
613     CALL write_ncdf_var('glamf',dimnames(1:2),name,Grid%glamf,'double')
614     CALL write_ncdf_var('gphit',dimnames(1:2),name,Grid%gphit,'double')
615     CALL write_ncdf_var('gphiu',dimnames(1:2),name,Grid%gphiu,'double')
616     CALL write_ncdf_var('gphiv',dimnames(1:2),name,Grid%gphiv,'double')
617     CALL write_ncdf_var('gphif',dimnames(1:2),name,Grid%gphif,'double')     
618
619     CALL write_ncdf_var('e1t',dimnames(1:2),name,Grid%e1t,'double')     
620     CALL write_ncdf_var('e1u',dimnames(1:2),name,Grid%e1u,'double')     
621     CALL write_ncdf_var('e1v',dimnames(1:2),name,Grid%e1v,'double')     
622     CALL write_ncdf_var('e1f',dimnames(1:2),name,Grid%e1f,'double')
623     CALL write_ncdf_var('e2t',dimnames(1:2),name,Grid%e2t,'double')
624     CALL write_ncdf_var('e2u',dimnames(1:2),name,Grid%e2u,'double')
625     CALL write_ncdf_var('e2v',dimnames(1:2),name,Grid%e2v,'double')
626     CALL write_ncdf_var('e2f',dimnames(1:2),name,Grid%e2f,'double')
627
628     CALL write_ncdf_var('ff_f',dimnames(1:2),name,ff_f,'double')
629     CALL write_ncdf_var('ff_t',dimnames(1:2),name,ff_t,'double')
630
631     CALL write_ncdf_var('e3t_1d',dimnames(3),name,e3t_1d,'double')
632     CALL write_ncdf_var('e3w_1d',dimnames(3),name,e3w_1d,'double')
633
634     CALL write_ncdf_var('e3t_0' ,dimnames(:),name,e3t_0 ,'double')
635     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double')
636     CALL write_ncdf_var('e3u_0' ,dimnames(:),name,e3u_0 ,'double')
637     CALL write_ncdf_var('e3v_0' ,dimnames(:),name,e3v_0 ,'double')
638     CALL write_ncdf_var('e3f_0' ,dimnames(:),name,e3f_0 ,'double')
639     CALL write_ncdf_var('e3w_0' ,dimnames(:),name,e3w_0 ,'double')
640     CALL write_ncdf_var('e3uw_0',dimnames(:),name,e3uw_0,'double')
641     CALL write_ncdf_var('e3vw_0',dimnames(:),name,e3vw_0,'double')
642
643     CALL write_ncdf_var('bottom_level',dimnames(1:2),name,bottom_level,'integer')
644     CALL write_ncdf_var('top_level'   ,dimnames(1:2),name,top_level   ,'integer')
645     CALL write_ncdf_var('bathy_meter' ,dimnames(1:2),name,zbathy      ,'float')
646
647     ! some attributes     
648     CALL copy_ncdf_att('nav_lon',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lon),MAXVAL(Grid%nav_lon))
649     CALL copy_ncdf_att('nav_lat',TRIM(parent_coordinate_file),name,MINVAL(Grid%nav_lat),MAXVAL(Grid%nav_lat))
650
651     CALL copy_ncdf_att('glamt',TRIM(parent_coordinate_file),name)
652     CALL copy_ncdf_att('glamu',TRIM(parent_coordinate_file),name)
653     CALL copy_ncdf_att('glamv',TRIM(parent_coordinate_file),name)
654     CALL copy_ncdf_att('glamf',TRIM(parent_coordinate_file),name)
655     CALL copy_ncdf_att('gphit',TRIM(parent_coordinate_file),name)
656     CALL copy_ncdf_att('gphiu',TRIM(parent_coordinate_file),name)
657     CALL copy_ncdf_att('gphiv',TRIM(parent_coordinate_file),name)
658     CALL copy_ncdf_att('gphif',TRIM(parent_coordinate_file),name)
659
660     CALL copy_ncdf_att('e1t',TRIM(parent_coordinate_file),name)
661     CALL copy_ncdf_att('e1u',TRIM(parent_coordinate_file),name)
662     CALL copy_ncdf_att('e1v',TRIM(parent_coordinate_file),name)
663     CALL copy_ncdf_att('e1f',TRIM(parent_coordinate_file),name)
664     CALL copy_ncdf_att('e2t',TRIM(parent_coordinate_file),name)
665     CALL copy_ncdf_att('e2u',TRIM(parent_coordinate_file),name)
666     CALL copy_ncdf_att('e2v',TRIM(parent_coordinate_file),name)
667     CALL copy_ncdf_att('e2f',TRIM(parent_coordinate_file),name)           
668     !
669     ! control print
670     WRITE(*,*) ' '
671     WRITE(*,*) 'Writing domcfg file: ',name
672     WRITE(*,*) ' '
673     !
674     DEALLOCATE( ff_t, ff_f, zbathy )
675     DEALLOCATE( bottom_level, top_level, mbathy )
676     DEALLOCATE( e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )
677     !     
678     write_domcfg = 1
679
680  END FUNCTION write_domcfg
681  !
682  SUBROUTINE zgr_z(e3t_1d, e3w_1d, gdept_1d, gdepw_1d )
683     !!----------------------------------------------------------------------
684     !!                   ***  ROUTINE zgr_z (from NEMO4) ***
685     !!                   
686     !! ** Purpose :   set the depth of model levels and the resulting
687     !!      vertical scale factors.
688     !!
689     !! ** Method  :   z-coordinate system (use in all type of coordinate)
690     !!        The depth of model levels is defined from an analytical
691     !!      function the derivative of which gives the scale factors.
692     !!        both depth and scale factors only depend on k (1d arrays).
693     !!              w-level: gdepw_1d  = gdep(k)
694     !!                       e3w_1d(k) = dk(gdep)(k)     = e3(k)
695     !!              t-level: gdept_1d  = gdep(k+0.5)
696     !!                       e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
697     !!
698     !! ** Action  : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
699     !!              - e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m)
700     !!
701     !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
702     !!----------------------------------------------------------------------
703     INTEGER  ::   jk                   ! dummy loop indices
704     REAL*8 ::   zt, zw                 ! temporary scalars
705     REAL*8 ::   zsur, za0, za1, zkth   ! Values set from parameters in
706     REAL*8 ::   zacr, zdzmin, zhmax    ! par_CONFIG_Rxx.h90
707     REAL*8 ::   za2, zkth2, zacr2      ! Values for optional double tanh function set from parameters
708
709     REAL*8, DIMENSION(:), INTENT(out) ::  e3t_1d, e3w_1d, gdept_1d, gdepw_1d
710     !!----------------------------------------------------------------------
711     !
712     !
713     ! Set variables from parameters
714     ! ------------------------------
715     zkth = ppkth       ;   zacr = ppacr
716     zdzmin = ppdzmin   ;   zhmax = pphmax
717     zkth2 = ppkth2     ;   zacr2 = ppacr2   ! optional (ldbletanh=T) double tanh parameters
718
719     ! If pa1 and pa0 and psur are et to pp_to_be_computed
720     !  za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
721    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
722         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN
723        !
724        za1  = (  ppdzmin - pphmax / FLOAT(N-1)  )                                                      &
725           & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(N-1) * (  LOG( COSH( (N - ppkth) / ppacr) )      &
726           &                                                   - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
727        za0  = ppdzmin - za1 *              TANH( (1-ppkth) / ppacr )
728        zsur =   - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
729     ELSE
730        za1 = pa1 ;       za0 = pa0 ;          zsur = psur
731        za2 = pa2                            ! optional (ldbletanh=T) double tanh parameter
732     ENDIF
733
734     ! Reference z-coordinate (depth - scale factor at T- and W-points)
735     ! ======================
736     IF( ppkth == 0. ) THEN            !  uniform vertical grid
737
738        za1 = zhmax / FLOAT(N-1)
739
740        DO jk = 1, N
741           zw = FLOAT( jk )
742           zt = FLOAT( jk ) + 0.5
743           gdepw_1d(jk) = ( zw - 1 ) * za1
744           gdept_1d(jk) = ( zt - 1 ) * za1
745           e3w_1d  (jk) =  za1
746           e3t_1d  (jk) =  za1
747        END DO
748     ELSE                                ! Madec & Imbard 1996 function
749        IF( .NOT. ldbletanh ) THEN
750           DO jk = 1, N
751              zw = REAL( jk )
752              zt = REAL( jk ) + 0.5
753              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) )  )
754              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) )  )
755              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth) / zacr   )
756              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth) / zacr   )
757           END DO
758        ELSE
759           DO jk = 1, N
760              zw = FLOAT( jk )
761              zt = FLOAT( jk ) + 0.5
762              ! Double tanh function
763              gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr  ) )    &
764                 &                             + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) )  )
765              gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr  ) )    &
766                 &                             + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) )  )
767              e3w_1d  (jk) =          za0      + za1        * TANH(       (zw-zkth ) / zacr  )      &
768                 &                             + za2        * TANH(       (zw-zkth2) / zacr2 )
769              e3t_1d  (jk) =          za0      + za1        * TANH(       (zt-zkth ) / zacr  )      &
770                 &                             + za2        * TANH(       (zt-zkth2) / zacr2 )
771           END DO
772        ENDIF
773        gdepw_1d(1) = 0.                    ! force first w-level to be exactly at zero
774     ENDIF
775
776     IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
777        !
778        !==>>>   need to be like this to compute the pressure gradient with ISF.
779        !        If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
780        !        define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
781        !
782        DO jk = 1, N-1
783           e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
784        END DO
785        e3t_1d(N) = e3t_1d(N-1)   ! we don't care because this level is masked in NEMO
786
787        DO jk = 2, N
788           e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
789        END DO
790        e3w_1d(1  ) = 2. * (gdept_1d(1) - gdepw_1d(1))
791     END IF
792
793     !
794  END SUBROUTINE zgr_z
795
796  SUBROUTINE zgr_zps( nx, ny, gdept_1d, gdepw_1d, e3t_1d, e3w_1d, zbathy, &
797     &                mbathy, e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0 )
798      !!----------------------------------------------------------------------
799      !!                  ***  ROUTINE zgr_zps (from NEMO4)  ***
800      !!                     
801      !! ** Purpose :   the depth and vertical scale factor in partial step
802      !!              reference z-coordinate case
803      !!
804      !! ** Method  :   Partial steps : computes the 3D vertical scale factors
805      !!      of T-, U-, V-, W-, UW-, VW and F-points that are associated with
806      !!      a partial step representation of bottom topography.
807      !!
808      !!        The reference depth of model levels is defined from an analytical
809      !!      function the derivative of which gives the reference vertical
810      !!      scale factors.
811      !!        From  depth and scale factors reference, we compute there new value
812      !!      with partial steps  on 3d arrays ( i, j, k ).
813      !!
814      !!              w-level: gdepw_0(i,j,k)  = gdep(k)
815      !!                       e3w_0(i,j,k) = dk(gdep)(k)     = e3(i,j,k)
816      !!              t-level: gdept_0(i,j,k)  = gdep(k+0.5)
817      !!                       e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
818      !!
819      !!        With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
820      !!      we find the mbathy index of the depth at each grid point.
821      !!      This leads us to three cases:
822      !!
823      !!              - bathy = 0 => mbathy = 0
824      !!              - 1 < mbathy < jpkm1   
825      !!              - bathy > gdepw_0(jpk) => mbathy = jpkm1 
826      !!
827      !!        Then, for each case, we find the new depth at t- and w- levels
828      !!      and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
829      !!      and f-points.
830      !!
831      !!        This routine is given as an example, it must be modified
832      !!      following the user s desiderata. nevertheless, the output as
833      !!      well as the way to compute the model levels and scale factors
834      !!      must be respected in order to insure second order accuracy
835      !!      schemes.
836      !!
837      !!         c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
838      !!         - - - - - - -   gdept_0, gdepw_0 and e3. are positives
839      !!     
840      !!  Reference :   Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
841      !!----------------------------------------------------------------------
842      INTEGER  ::   ji, jj, jk       ! dummy loop indices
843      INTEGER  ::   ik
844      INTEGER  ::   ikb, ikt       ! temporary integers
845      REAL*8 ::   ze3tp , ze3wp    ! Last ocean level thickness at T- and W-points
846      REAL*8 ::   zdepwp, zdepth   ! Ajusted ocean depth to avoid too small e3t
847      REAL*8 ::   zmax             ! temporary scalar
848      !
849      REAL*8, ALLOCATABLE, DIMENSION(:,:,:) ::  gdept_0, gdepw_0, zprt
850      !
851      INTEGER,                   INTENT(in   ) ::  nx, ny
852      REAL*8 , DIMENSION(:)    , INTENT(in   ) ::  gdept_1d, gdepw_1d, e3t_1d, e3w_1d
853      REAL*8 , DIMENSION(:,:)  , INTENT(inout) ::  zbathy
854      INTEGER, DIMENSION(:,:)  , INTENT(  out) ::  mbathy
855      REAL*8 , DIMENSION(:,:,:), INTENT(  out) ::  e3t_0, e3u_0, e3v_0, e3f_0, e3w_0, e3uw_0, e3vw_0
856      !
857      !!---------------------------------------------------------------------
858     
859      ALLOCATE( zprt(nx,ny,N), gdept_0(nx,ny,N), gdepw_0(nx,ny,N) )
860      !
861      ! bathymetry in level (from bathy_meter)
862      ! ===================
863      zmax = gdepw_1d(N) + e3t_1d(N)        ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(N-1) )
864      zbathy(:,:) = MIN( zmax ,  zbathy(:,:) )    ! bounded value of bathy (min already set at the end of zgr_bat)
865      WHERE( zbathy(:,:) == 0. )   ;   mbathy(:,:) = 0       ! land  : set mbathy to 0
866      ELSEWHERE                    ;   mbathy(:,:) = N-1   ! ocean : initialize mbathy to the max ocean level
867      END WHERE
868
869      ! Compute mbathy for ocean points (i.e. the number of ocean levels)
870      ! find the number of ocean levels such that the last level thickness
871      ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
872      ! e3t_1d is the reference level thickness
873      DO jk = N-1, 1, -1
874         zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
875         WHERE( 0. < zbathy(:,:) .AND. zbathy(:,:) <= zdepth )   mbathy(:,:) = jk-1
876      END DO
877
878      ! Scale factors and depth at T- and W-points
879      DO jk = 1, N                        ! intitialization to the reference z-coordinate
880         gdept_0(:,:,jk) = gdept_1d(jk)
881         gdepw_0(:,:,jk) = gdepw_1d(jk)
882         e3t_0  (:,:,jk) = e3t_1d  (jk)
883         e3w_0  (:,:,jk) = e3w_1d  (jk)
884      END DO
885     
886      ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
887      !!clem: not implemented yet
888      !!      IF ( ln_isfcav == 1 ) CALL zgr_isf
889      !
890      ! Scale factors and depth at T- and W-points
891!      IF ( ln_isfcav == 0 ) THEN
892         DO jj = 1, ny
893            DO ji = 1, nx
894               ik = mbathy(ji,jj)
895               IF( ik > 0 ) THEN               ! ocean point only
896                  ! max ocean level case
897                  IF( ik == N-1 ) THEN
898                     zdepwp = zbathy(ji,jj)
899                     ze3tp  = zbathy(ji,jj) - gdepw_1d(ik)
900                     ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) )
901                     e3t_0(ji,jj,ik  ) = ze3tp
902                     e3t_0(ji,jj,ik+1) = ze3tp
903                     e3w_0(ji,jj,ik  ) = ze3wp
904                     e3w_0(ji,jj,ik+1) = ze3tp
905                     gdepw_0(ji,jj,ik+1) = zdepwp
906                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
907                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
908                     !
909                  ELSE                         ! standard case
910                     IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = zbathy(ji,jj)
911                     ELSE                                        ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
912                     ENDIF
913   !gm Bug?  check the gdepw_1d
914                     !       ... on ik
915                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
916                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
917                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
918                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
919                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
920                     e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) )   &
921                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
922                     !       ... on ik+1
923                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
924                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
925                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
926                  ENDIF
927               ENDIF
928            END DO
929         END DO
930         !
931!         DO jj = 1, ny
932!            DO ji = 1, nx
933!               ik = mbathy(ji,jj)
934!               IF( ik > 0 ) THEN               ! ocean point only
935!                  e3tp (ji,jj) = e3t_0(ji,jj,ik)
936!                  e3wp (ji,jj) = e3w_0(ji,jj,ik)
937!               ENDIF
938!            END DO
939!         END DO
940!      END IF
941      !
942      ! Scale factors and depth at U-, V-, UW and VW-points
943      DO jk = 1, N                        ! initialisation to z-scale factors
944         e3u_0 (:,:,jk) = e3t_1d(jk)
945         e3v_0 (:,:,jk) = e3t_1d(jk)
946         e3uw_0(:,:,jk) = e3w_1d(jk)
947         e3vw_0(:,:,jk) = e3w_1d(jk)
948      END DO
949
950      DO jk = 1,N                         ! Computed as the minimum of neighbooring scale factors
951         DO jj = 1, ny-1
952            DO ji = 1, nx-1   ! vector opt.
953               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
954               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
955               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
956               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
957            END DO
958         END DO
959      END DO
960     
961!      IF ( ln_isfcav == 1 ) THEN
962!      ! (ISF) define e3uw (adapted for 2 cells in the water column)
963!         DO jj = 2, ny-1
964!            DO ji = 2, nx-1   ! vector opt.
965!               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
966!               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
967!               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
968!                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
969!               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
970!               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
971!               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
972!                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
973!            END DO
974!         END DO
975!      END IF
976      !
977
978      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries)
979         WHERE( e3u_0 (:,:,jk) == 0. )   e3u_0 (:,:,jk) = e3t_1d(jk)
980         WHERE( e3v_0 (:,:,jk) == 0. )   e3v_0 (:,:,jk) = e3t_1d(jk)
981         WHERE( e3uw_0(:,:,jk) == 0. )   e3uw_0(:,:,jk) = e3w_1d(jk)
982         WHERE( e3vw_0(:,:,jk) == 0. )   e3vw_0(:,:,jk) = e3w_1d(jk)
983      END DO
984     
985      ! Scale factor at F-point
986      DO jk = 1, N                        ! initialisation to z-scale factors
987         e3f_0(:,:,jk) = e3t_1d(jk)
988      END DO
989      DO jk = 1, N                        ! Computed as the minimum of neighbooring V-scale factors
990         DO jj = 1, ny-1
991            DO ji = 1, nx-1   ! vector opt.
992               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
993            END DO
994         END DO
995      END DO
996      !
997      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries)
998         WHERE( e3f_0(:,:,jk) == 0. )   e3f_0(:,:,jk) = e3t_1d(jk)
999      END DO
1000
1001      ! Control of the sign
1002      IF( MINVAL( e3t_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3t_0 <= 0' 
1003      IF( MINVAL( e3w_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3w_0 <= 0' 
1004      IF( MINVAL( gdept_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdept_0 <  0' 
1005      IF( MINVAL( gdepw_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdepw_0 <  0' 
1006     
1007      ! Compute gde3w_0 (vertical sum of e3w)
1008!      IF ( ln_isfcav ==1 ) THEN ! if cavity
1009!         WHERE( misfdep == 0 )   misfdep = 1
1010!         DO jj = 1,ny
1011!            DO ji = 1,nx
1012!               gde3w_0(ji,jj,1) = 0.5 * e3w_0(ji,jj,1)
1013!               DO jk = 2, misfdep(ji,jj)
1014!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1015!               END DO
1016!               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5 * e3w_0(ji,jj,misfdep(ji,jj))
1017!               DO jk = misfdep(ji,jj) + 1, N
1018!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1019!               END DO
1020!            END DO
1021!         END DO
1022!      ELSE ! no cavity
1023!         gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
1024!         DO jk = 2, N
1025!            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1026!         END DO
1027!      END IF
1028      !
1029      DEALLOCATE( zprt, gdept_0, gdepw_0 )
1030      !
1031   END SUBROUTINE zgr_zps
1032
1033  !*****************************************************
1034  !   function set_child_name(Parentname,Childname)
1035  !*****************************************************
1036  !
1037  SUBROUTINE set_child_name(Parentname,Childname)
1038    !
1039    CHARACTER(*),INTENT(in) :: Parentname
1040    CHARACTER(*),INTENT(out) :: Childname
1041    CHARACTER(2) :: prefix
1042    INTEGER :: pos
1043    !   
1044    pos  = INDEX(TRIM(Parentname),'/',back=.TRUE.)
1045    !
1046    prefix=Parentname(pos+1:pos+2)
1047    IF (prefix == '1_') THEN
1048       Childname = '2_'//Parentname(pos+3:LEN(Parentname)) 
1049    ELSEIF (prefix == '2_') THEN
1050       Childname = '3_'//Parentname(pos+3:LEN(Parentname)) 
1051    ELSEIF (prefix == '3_') THEN
1052       Childname = '4_'//Parentname(pos+3:LEN(Parentname)) 
1053    ELSEIF (prefix == '4_') THEN
1054       Childname = '5_'//Parentname(pos+3:LEN(Parentname)) 
1055    ELSE
1056       Childname = '1_'//Parentname(pos+1:LEN(Parentname)) 
1057    ENDIF
1058    !   
1059  END SUBROUTINE set_child_name
1060  !
1061  !*****************************************************
1062  !   subroutine get_interptype(varname,interp_type,conservation)
1063  !*****************************************************
1064  !
1065  SUBROUTINE get_interptype( varname,interp_type,conservation )
1066    !
1067    LOGICAL,OPTIONAL :: conservation
1068    CHARACTER(*) :: interp_type,varname
1069    INTEGER :: pos,pos1,pos2,pos3,i,k 
1070    LOGICAL :: find       
1071    i=1
1072    DO WHILE ( TRIM(VAR_INTERP(i)) /= '' )     
1073       pos = INDEX( TRIM(VAR_INTERP(i)) , TRIM(varname) )
1074       IF ( pos .NE. 0 ) THEN     
1075          pos1 = INDEX( TRIM(VAR_INTERP(i)) , 'bicubic' )
1076          pos2 = INDEX( TRIM(VAR_INTERP(i)) , 'bilinear' )
1077          pos3 = INDEX( TRIM(VAR_INTERP(i)) , 'conservative' )
1078          ! initialize interp_type
1079          IF( pos1 .NE. 0 ) interp_type = 'bicubic'
1080          IF( pos2 .NE. 0 ) interp_type = 'bilinear'
1081          IF( pos1 .EQ. 0 .AND. pos2 .EQ. 0) interp_type = 'bicubic'
1082          ! initialize conservation                                           
1083          IF( pos3 .NE. 0 .AND. PRESENT(conservation) ) THEN
1084             conservation = .TRUE.
1085             RETURN
1086          ELSE
1087             conservation = .FALSE.
1088          ENDIF
1089          find = .FALSE.
1090          IF( PRESENT(conservation) ) THEN
1091             k=0
1092             conservation = .FALSE.         
1093             DO WHILE( k < SIZE(flxtab) .AND. .NOT.find ) 
1094                k = k+1
1095                IF( TRIM(varname) .EQ. TRIM(flxtab(k)) ) THEN       
1096                   conservation = .TRUE.
1097                   find = .TRUE.
1098                ENDIF
1099             END DO
1100          ENDIF
1101          RETURN
1102       ENDIF
1103       i = i+1
1104    END DO
1105    !default values interp_type = bicubic // conservation = false     
1106    interp_type = 'bicubic'     
1107    IF( PRESENT(conservation) ) conservation = .FALSE.
1108
1109    RETURN
1110    !   
1111  END SUBROUTINE get_interptype
1112  !     
1113  !*****************************************************
1114  !   subroutine Init_mask(name,Grid)
1115  !*****************************************************
1116  !
1117  SUBROUTINE Init_mask(name,Grid,jpiglo,jpjglo)
1118    !
1119    USE io_netcdf
1120    !     
1121    CHARACTER(*) name
1122    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
1123    TYPE(Coordinates) :: Grid
1124    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
1125    !
1126    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
1127       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
1128    ELSE
1129       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
1130    ENDIF
1131
1132
1133    !
1134    WRITE(*,*) 'Init masks in T,U,V,F points'   
1135    !
1136    nx = SIZE(Grid%Bathy_level,1)
1137    ny = SIZE(Grid%Bathy_level,2)
1138    !
1139    !     
1140    ALLOCATE(Grid%tmask(nx,ny,N), &
1141         Grid%umask(nx,ny,N), &
1142         Grid%vmask(nx,ny,N), &
1143         Grid%fmask(nx,ny,N))
1144    !
1145    DO k = 1,N
1146       !     
1147       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
1148          Grid%tmask(:,:,k) = 0
1149       ELSEWHERE
1150          Grid%tmask(:,:,k) = 1
1151       END WHERE
1152       !
1153    END DO
1154    !
1155    Grid%umask(1:nx-1,:,:) = Grid%tmask(1:nx-1,:,:)*Grid%tmask(2:nx,:,:)
1156    Grid%vmask(:,1:ny-1,:) = Grid%tmask(:,1:ny-1,:)*Grid%tmask(:,2:ny,:)
1157    !
1158    Grid%umask(nx,:,:) = Grid%tmask(nx,:,:)
1159    Grid%vmask(:,ny,:) = Grid%tmask(:,ny,:)
1160    !     
1161    Grid%fmask(1:nx-1,1:ny-1,:) = Grid%tmask(1:nx-1,1:ny-1,:)*Grid%tmask(2:nx,1:ny-1,:)* &
1162         Grid%tmask(1:nx-1,2:ny,:)*Grid%tmask(2:nx,2:ny,:) 
1163    !
1164    Grid%fmask(nx,:,:) = Grid%tmask(nx,:,:)
1165    Grid%fmask(:,ny,:) = Grid%tmask(:,ny,:)
1166    !
1167    ALLOCATE(zwf(nx,ny))
1168    !     
1169    DO k = 1,N
1170       !
1171       zwf(:,:) = Grid%fmask(:,:,k)
1172       !         
1173       DO j = 2, ny-1
1174          DO i = 2,nx-1           
1175             IF( Grid%fmask(i,j,k) == 0. ) THEN               
1176                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)))
1177             END IF
1178          END DO
1179       END DO
1180       !
1181       DO j = 2, ny-1
1182          IF( Grid%fmask(1,j,k) == 0. ) THEN
1183             Grid%fmask(1,j,k) = shlat * MIN(1.,MAX(zwf(2,j),zwf(1,j+1),zwf(1,j-1)))
1184          ENDIF
1185
1186          IF( Grid%fmask(nx,j,k) == 0. ) THEN
1187             Grid%fmask(nx,j,k) = shlat * MIN(1.,MAX(zwf(nx,j+1),zwf(nx-1,j),zwf(nx,j-1)))
1188          ENDIF
1189       END DO
1190       !         
1191       DO i = 2, nx-1         
1192          IF( Grid%fmask(i,1,k) == 0. ) THEN
1193             Grid%fmask(i, 1 ,k) = shlat*MIN( 1.,MAX(zwf(i+1,1),zwf(i,2),zwf(i-1,1)))
1194          ENDIF
1195          !           
1196          IF( Grid%fmask(i,ny,k) == 0. ) THEN
1197             Grid%fmask(i,ny,k) = shlat * MIN(1.,MAX(zwf(i+1,ny),zwf(i-1,ny),zwf(i,ny-1)))
1198          ENDIF
1199       END DO
1200       !!
1201    END DO
1202    !!     
1203  END SUBROUTINE Init_mask
1204  !
1205  !*****************************************************
1206  !   subroutine Init_Tmask(name,Grid)
1207  !*****************************************************
1208  !
1209  SUBROUTINE Init_Tmask(name,Grid,jpiglo,jpjglo)
1210    !
1211    USE io_netcdf
1212    !     
1213    CHARACTER(*) name
1214    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
1215    TYPE(Coordinates) :: Grid
1216    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
1217    !
1218    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
1219       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
1220    ELSE
1221       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
1222    ENDIF
1223    !
1224    nx = SIZE(Grid%Bathy_level,1)
1225    ny = SIZE(Grid%Bathy_level,2) 
1226    !
1227    WRITE(*,*) 'Init masks in T points'   
1228    !     
1229    ALLOCATE( Grid%tmask(nx,ny,N) )
1230    !
1231    DO k = 1,N
1232       !     
1233       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
1234          Grid%tmask(:,:,k) = 0.
1235       ELSEWHERE
1236          Grid%tmask(:,:,k) = 1.
1237       END WHERE
1238       !
1239    END DO
1240    !     
1241  END SUBROUTINE Init_Tmask
1242  !
1243  !*****************************************************
1244  !   subroutine Init_ssmask(name,Grid)
1245  !*****************************************************
1246  !
1247!  SUBROUTINE Init_ssmask(varname,filename,Grid,jpiglo,jpjglo)
1248!    !
1249!    USE io_netcdf
1250!    !     
1251!    CHARACTER(*) varname,filename
1252!    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
1253!    TYPE(Coordinates) :: Grid
1254!    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
1255!    !
1256!    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
1257!       CALL read_ncdf_var(varname,filename,Grid%bathy_level)
1258!    ELSE
1259!       CALL read_ncdf_var(varname,filename,Grid%bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
1260!    ENDIF
1261!    !
1262!    nx = SIZE(Grid%bathy_level,1)
1263!    ny = SIZE(Grid%bathy_level,2)
1264!    !
1265!    WRITE(*,*) 'Init surface masks in T points'   
1266!    !     
1267!    ALLOCATE( Grid%ssmask(nx,ny), Grid%ssumask(nx,ny), Grid%ssvmask(nx,ny) )
1268!    !
1269!    WHERE(Grid%bathy_level(:,:) <= 0. )   
1270!       Grid%ssmask(:,:) = 0.
1271!    ELSEWHERE
1272!       Grid%ssmask(:,:) = 1.
1273!    END WHERE
1274!    !
1275!    Grid%ssumask(1:nx-1,:) = Grid%ssmask(1:nx-1,:)*Grid%ssmask(2:nx,:)
1276!    Grid%ssvmask(:,1:ny-1) = Grid%ssmask(:,1:ny-1)*Grid%ssmask(:,2:ny)
1277!    !
1278!    Grid%ssumask(nx,:) = Grid%ssmask(nx,:)
1279!    Grid%ssvmask(:,ny) = Grid%ssmask(:,ny)
1280!    !     
1281!  END SUBROUTINE Init_ssmask
1282  !
1283  !*****************************************************
1284  !   subroutine get_mask(name,Grid)
1285  !*****************************************************
1286  !
1287  SUBROUTINE get_mask(level,posvar,mask,filename)
1288    !
1289    USE io_netcdf
1290    !     
1291    CHARACTER(*) filename
1292    CHARACTER(*) posvar
1293    INTEGER :: level, nx, ny
1294    LOGICAL,DIMENSION(:,:),POINTER :: mask
1295    INTEGER,DIMENSION(:,:),POINTER :: maskT,maskU,maskV
1296    !     
1297    TYPE(Coordinates) :: Grid
1298    !
1299    CALL read_ncdf_var('Bathy_level',filename,Grid%Bathy_level)
1300    !
1301    nx = SIZE(Grid%Bathy_level,1)
1302    ny = SIZE(Grid%Bathy_level,2)
1303    ALLOCATE(maskT(nx,ny),mask(nx,ny))
1304    mask = .TRUE.
1305    !
1306    WHERE(Grid%Bathy_level(:,:) <= level-1 )   
1307       maskT(:,:) = 0
1308    ELSEWHERE
1309       maskT(:,:) = 1
1310    END WHERE
1311    !
1312    SELECT CASE(posvar)
1313       !
1314    CASE('T')
1315       !
1316       WHERE(maskT > 0)
1317          mask = .TRUE.
1318       ELSEWHERE
1319          mask = .FALSE.
1320       END WHERE
1321       DEALLOCATE(maskT)
1322       !
1323    CASE('U')
1324       !
1325       ALLOCATE(maskU(nx,ny))
1326       maskU(1:nx-1,:) = maskT(1:nx-1,:)*maskT(2:nx,:)
1327       maskU(nx,:) = maskT(nx,:)
1328       WHERE(maskU > 0)
1329          mask = .TRUE.
1330       ELSEWHERE
1331          mask = .FALSE.
1332       END WHERE
1333       DEALLOCATE(maskU,maskT)
1334       !
1335    CASE('V')
1336       !
1337       ALLOCATE(maskV(nx,ny))
1338       maskV(:,1:ny-1) = maskT(:,1:ny-1)*maskT(:,2:ny)
1339       maskV(:,ny) = maskT(:,ny)
1340       WHERE(maskT > 0)
1341          mask = .TRUE.
1342       ELSEWHERE
1343          mask = .FALSE.
1344       END WHERE
1345       DEALLOCATE(maskV,maskT)
1346       !
1347    END SELECT
1348    !     
1349  END SUBROUTINE get_mask
1350  !
1351  !
1352  !*****************************************************
1353  !   subroutine read_dimg_var(unit,irec,field)
1354  !*****************************************************
1355  !
1356  SUBROUTINE read_dimg_var(unit,irec,field,jpk)
1357    !     
1358    INTEGER :: unit,irec,jpk
1359    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
1360    INTEGER :: k
1361    !     
1362    DO k = 1,jpk
1363       READ(unit,REC=irec) field(:,:,k,1)
1364       irec = irec + 1
1365    ENDDO
1366    !
1367  END SUBROUTINE read_dimg_var
1368  !
1369  !
1370  !*****************************************************
1371  !   subroutine read_dimg_var(unit,irec,field)
1372  !*****************************************************
1373  !
1374  SUBROUTINE write_dimg_var(unit,irec,field,jpk)
1375    !     
1376    INTEGER :: unit,irec,jpk
1377    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
1378    INTEGER :: k
1379    !     
1380    DO k = 1,jpk
1381       WRITE(unit,REC=irec) field(:,:,k,1)
1382       irec = irec + 1
1383    ENDDO
1384    !
1385  END SUBROUTINE write_dimg_var
1386
1387END MODULE agrif_readwrite
Note: See TracBrowser for help on using the repository browser.