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

Last change on this file was 12264, checked in by clem, 14 months ago

finishing cleaning the nesting tools (as far as I could). I will stop there

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