New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_readwrite.f90 in utils/tools/NESTING/src – NEMO

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

Last change on this file since 12253 was 12253, checked in by clem, 4 years ago

1) resolve some conflicts when using partial steps or not. 2) Make sure grids remain identical when they should. 3) in case of a double zoom, it is no more necessary to have the bilinear interpolation option when updating the 1st parent grid. Btw, I know these comments are unclear but it is very difficult to explain...
  • Property svn:keywords set to Id
File size: 53.4 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      IF( partial_steps ) THEN
900         ! Scale factors and depth at T- and W-points
901         ! IF ( ln_isfcav == 0 ) THEN
902         DO jj = 1, ny
903            DO ji = 1, nx
904               ik = mbathy(ji,jj)
905               IF( ik > 0 ) THEN               ! ocean point only
906                  ! max ocean level case
907                  IF( ik == N-1 ) THEN
908                     zdepwp = zbathy(ji,jj)
909                     ze3tp  = zbathy(ji,jj) - gdepw_1d(ik)
910                     ze3wp = 0.5 * e3w_1d(ik) * ( 1. + ( ze3tp/e3t_1d(ik) ) )
911                     e3t_0(ji,jj,ik  ) = ze3tp
912                     e3t_0(ji,jj,ik+1) = ze3tp
913                     e3w_0(ji,jj,ik  ) = ze3wp
914                     e3w_0(ji,jj,ik+1) = ze3tp
915                     gdepw_0(ji,jj,ik+1) = zdepwp
916                     gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp
917                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
918                     !
919                  ELSE                         ! standard case
920                     IF( zbathy(ji,jj) <= gdepw_1d(ik+1) ) THEN  ;   gdepw_0(ji,jj,ik+1) = zbathy(ji,jj)
921                     ELSE                                        ;   gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
922                     ENDIF
923                     !gm Bug?  check the gdepw_1d
924                     !       ... on ik
925                     gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) )   &
926                        &                             * ((gdept_1d(     ik  ) - gdepw_1d(ik) )   &
927                        &                             / ( gdepw_1d(     ik+1) - gdepw_1d(ik) ))
928                     e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) )   & 
929                        &                             / ( gdepw_1d(      ik+1) - gdepw_1d(ik) ) 
930                     e3w_0(ji,jj,ik) = 0.5 * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2. * gdepw_1d(ik) )   &
931                        &                     * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
932                     !       ... on ik+1
933                     e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
934                     e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik)
935                     gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
936                  ENDIF
937               ENDIF
938            END DO
939         END DO
940         !
941         ! DO jj = 1, ny
942         !    DO ji = 1, nx
943         !       ik = mbathy(ji,jj)
944         !       IF( ik > 0 ) THEN               ! ocean point only
945         !          e3tp (ji,jj) = e3t_0(ji,jj,ik)
946         !          e3wp (ji,jj) = e3w_0(ji,jj,ik)
947         !       ENDIF
948         !    END DO
949         ! END DO
950         ! END IF
951      ENDIF
952      !
953      ! Scale factors and depth at U-, V-, UW and VW-points
954      DO jk = 1, N                        ! initialisation to z-scale factors
955         e3u_0 (:,:,jk) = e3t_1d(jk)
956         e3v_0 (:,:,jk) = e3t_1d(jk)
957         e3uw_0(:,:,jk) = e3w_1d(jk)
958         e3vw_0(:,:,jk) = e3w_1d(jk)
959      END DO
960
961      DO jk = 1,N                         ! Computed as the minimum of neighbooring scale factors
962         DO jj = 1, ny-1
963            DO ji = 1, nx-1   ! vector opt.
964               e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
965               e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
966               e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
967               e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
968            END DO
969         END DO
970      END DO
971     
972!      IF ( ln_isfcav == 1 ) THEN
973!      ! (ISF) define e3uw (adapted for 2 cells in the water column)
974!         DO jj = 2, ny-1
975!            DO ji = 2, nx-1   ! vector opt.
976!               ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
977!               ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
978!               IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji+1,jj  ,ikb  ) ) &
979!                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj  ,ikb-1) )
980!               ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
981!               ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
982!               IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) =  MIN( gdept_0(ji,jj,ikb  ), gdept_0(ji  ,jj+1,ikb  ) ) &
983!                                       &            - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji  ,jj+1,ikb-1) )
984!            END DO
985!         END DO
986!      END IF
987      !
988
989      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries)
990         WHERE( e3u_0 (:,:,jk) == 0. )   e3u_0 (:,:,jk) = e3t_1d(jk)
991         WHERE( e3v_0 (:,:,jk) == 0. )   e3v_0 (:,:,jk) = e3t_1d(jk)
992         WHERE( e3uw_0(:,:,jk) == 0. )   e3uw_0(:,:,jk) = e3w_1d(jk)
993         WHERE( e3vw_0(:,:,jk) == 0. )   e3vw_0(:,:,jk) = e3w_1d(jk)
994      END DO
995     
996      ! Scale factor at F-point
997      DO jk = 1, N                        ! initialisation to z-scale factors
998         e3f_0(:,:,jk) = e3t_1d(jk)
999      END DO
1000      DO jk = 1, N                        ! Computed as the minimum of neighbooring V-scale factors
1001         DO jj = 1, ny-1
1002            DO ji = 1, nx-1   ! vector opt.
1003               e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1004            END DO
1005         END DO
1006      END DO
1007      !
1008      DO jk = 1, N                        ! set to z-scale factor if zero (i.e. along closed boundaries)
1009         WHERE( e3f_0(:,:,jk) == 0. )   e3f_0(:,:,jk) = e3t_1d(jk)
1010      END DO
1011
1012      ! Control of the sign
1013      IF( MINVAL( e3t_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3t_0 <= 0' 
1014      IF( MINVAL( e3w_0  (:,:,:) ) <= 0. )   STOP '    zgr_zps :   e r r o r   e3w_0 <= 0' 
1015      IF( MINVAL( gdept_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdept_0 <  0' 
1016      IF( MINVAL( gdepw_0(:,:,:) ) <  0. )   STOP '    zgr_zps :   e r r o r   gdepw_0 <  0' 
1017     
1018      ! Compute gde3w_0 (vertical sum of e3w)
1019!      IF ( ln_isfcav ==1 ) THEN ! if cavity
1020!         WHERE( misfdep == 0 )   misfdep = 1
1021!         DO jj = 1,ny
1022!            DO ji = 1,nx
1023!               gde3w_0(ji,jj,1) = 0.5 * e3w_0(ji,jj,1)
1024!               DO jk = 2, misfdep(ji,jj)
1025!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1026!               END DO
1027!               IF( misfdep(ji,jj) >= 2 )   gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5 * e3w_0(ji,jj,misfdep(ji,jj))
1028!               DO jk = misfdep(ji,jj) + 1, N
1029!                  gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
1030!               END DO
1031!            END DO
1032!         END DO
1033!      ELSE ! no cavity
1034!         gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
1035!         DO jk = 2, N
1036!            gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1037!         END DO
1038!      END IF
1039      !
1040      DEALLOCATE( zprt, gdept_0, gdepw_0 )
1041      !
1042   END SUBROUTINE zgr_zps
1043
1044  !*****************************************************
1045  !   function set_child_name(Parentname,Childname)
1046  !*****************************************************
1047  !
1048  SUBROUTINE set_child_name(Parentname,Childname)
1049    !
1050    CHARACTER(*),INTENT(in) :: Parentname
1051    CHARACTER(*),INTENT(out) :: Childname
1052    CHARACTER(2) :: prefix
1053    INTEGER :: pos
1054    !   
1055    pos  = INDEX(TRIM(Parentname),'/',back=.TRUE.)
1056    !
1057    prefix=Parentname(pos+1:pos+2)
1058    IF (prefix == '1_') THEN
1059       Childname = '2_'//Parentname(pos+3:LEN(Parentname)) 
1060    ELSEIF (prefix == '2_') THEN
1061       Childname = '3_'//Parentname(pos+3:LEN(Parentname)) 
1062    ELSEIF (prefix == '3_') THEN
1063       Childname = '4_'//Parentname(pos+3:LEN(Parentname)) 
1064    ELSEIF (prefix == '4_') THEN
1065       Childname = '5_'//Parentname(pos+3:LEN(Parentname)) 
1066    ELSE
1067       Childname = '1_'//Parentname(pos+1:LEN(Parentname)) 
1068    ENDIF
1069    !   
1070  END SUBROUTINE set_child_name
1071  !
1072  !*****************************************************
1073  !   subroutine get_interptype(varname,interp_type,conservation)
1074  !*****************************************************
1075  !
1076  SUBROUTINE get_interptype( varname,interp_type,conservation )
1077    !
1078    LOGICAL,OPTIONAL :: conservation
1079    CHARACTER(*) :: interp_type,varname
1080    INTEGER :: pos,pos1,pos2,pos3,i,k 
1081    LOGICAL :: find       
1082    i=1
1083    DO WHILE ( TRIM(VAR_INTERP(i)) /= '' )     
1084       pos = INDEX( TRIM(VAR_INTERP(i)) , TRIM(varname) )
1085       IF ( pos .NE. 0 ) THEN     
1086          pos1 = INDEX( TRIM(VAR_INTERP(i)) , 'bicubic' )
1087          pos2 = INDEX( TRIM(VAR_INTERP(i)) , 'bilinear' )
1088          pos3 = INDEX( TRIM(VAR_INTERP(i)) , 'conservative' )
1089          ! initialize interp_type
1090          IF( pos1 .NE. 0 ) interp_type = 'bicubic'
1091          IF( pos2 .NE. 0 ) interp_type = 'bilinear'
1092          IF( pos1 .EQ. 0 .AND. pos2 .EQ. 0) interp_type = 'bicubic'
1093          ! initialize conservation                                           
1094          IF( pos3 .NE. 0 .AND. PRESENT(conservation) ) THEN
1095             conservation = .TRUE.
1096             RETURN
1097          ELSE
1098             conservation = .FALSE.
1099          ENDIF
1100          find = .FALSE.
1101          IF( PRESENT(conservation) ) THEN
1102             k=0
1103             conservation = .FALSE.         
1104             DO WHILE( k < SIZE(flxtab) .AND. .NOT.find ) 
1105                k = k+1
1106                IF( TRIM(varname) .EQ. TRIM(flxtab(k)) ) THEN       
1107                   conservation = .TRUE.
1108                   find = .TRUE.
1109                ENDIF
1110             END DO
1111          ENDIF
1112          RETURN
1113       ENDIF
1114       i = i+1
1115    END DO
1116    !default values interp_type = bicubic // conservation = false     
1117    interp_type = 'bicubic'     
1118    IF( PRESENT(conservation) ) conservation = .FALSE.
1119
1120    RETURN
1121    !   
1122  END SUBROUTINE get_interptype
1123  !     
1124  !*****************************************************
1125  !   subroutine Init_mask(name,Grid)
1126  !*****************************************************
1127  !
1128  SUBROUTINE Init_mask(name,Grid,jpiglo,jpjglo)
1129    !
1130    USE io_netcdf
1131    !     
1132    CHARACTER(*) name
1133    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
1134    TYPE(Coordinates) :: Grid
1135    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
1136    !
1137    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
1138       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
1139    ELSE
1140       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
1141    ENDIF
1142
1143
1144    !
1145    WRITE(*,*) 'Init masks in T,U,V,F points'   
1146    !
1147    nx = SIZE(Grid%Bathy_level,1)
1148    ny = SIZE(Grid%Bathy_level,2)
1149    !
1150    !     
1151    ALLOCATE(Grid%tmask(nx,ny,N), &
1152         Grid%umask(nx,ny,N), &
1153         Grid%vmask(nx,ny,N), &
1154         Grid%fmask(nx,ny,N))
1155    !
1156    DO k = 1,N
1157       !     
1158       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
1159          Grid%tmask(:,:,k) = 0
1160       ELSEWHERE
1161          Grid%tmask(:,:,k) = 1
1162       END WHERE
1163       !
1164    END DO
1165    !
1166    Grid%umask(1:nx-1,:,:) = Grid%tmask(1:nx-1,:,:)*Grid%tmask(2:nx,:,:)
1167    Grid%vmask(:,1:ny-1,:) = Grid%tmask(:,1:ny-1,:)*Grid%tmask(:,2:ny,:)
1168    !
1169    Grid%umask(nx,:,:) = Grid%tmask(nx,:,:)
1170    Grid%vmask(:,ny,:) = Grid%tmask(:,ny,:)
1171    !     
1172    Grid%fmask(1:nx-1,1:ny-1,:) = Grid%tmask(1:nx-1,1:ny-1,:)*Grid%tmask(2:nx,1:ny-1,:)* &
1173         Grid%tmask(1:nx-1,2:ny,:)*Grid%tmask(2:nx,2:ny,:) 
1174    !
1175    Grid%fmask(nx,:,:) = Grid%tmask(nx,:,:)
1176    Grid%fmask(:,ny,:) = Grid%tmask(:,ny,:)
1177    !
1178    ALLOCATE(zwf(nx,ny))
1179    !     
1180    DO k = 1,N
1181       !
1182       zwf(:,:) = Grid%fmask(:,:,k)
1183       !         
1184       DO j = 2, ny-1
1185          DO i = 2,nx-1           
1186             IF( Grid%fmask(i,j,k) == 0. ) THEN               
1187                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)))
1188             END IF
1189          END DO
1190       END DO
1191       !
1192       DO j = 2, ny-1
1193          IF( Grid%fmask(1,j,k) == 0. ) THEN
1194             Grid%fmask(1,j,k) = shlat * MIN(1.,MAX(zwf(2,j),zwf(1,j+1),zwf(1,j-1)))
1195          ENDIF
1196
1197          IF( Grid%fmask(nx,j,k) == 0. ) THEN
1198             Grid%fmask(nx,j,k) = shlat * MIN(1.,MAX(zwf(nx,j+1),zwf(nx-1,j),zwf(nx,j-1)))
1199          ENDIF
1200       END DO
1201       !         
1202       DO i = 2, nx-1         
1203          IF( Grid%fmask(i,1,k) == 0. ) THEN
1204             Grid%fmask(i, 1 ,k) = shlat*MIN( 1.,MAX(zwf(i+1,1),zwf(i,2),zwf(i-1,1)))
1205          ENDIF
1206          !           
1207          IF( Grid%fmask(i,ny,k) == 0. ) THEN
1208             Grid%fmask(i,ny,k) = shlat * MIN(1.,MAX(zwf(i+1,ny),zwf(i-1,ny),zwf(i,ny-1)))
1209          ENDIF
1210       END DO
1211       !!
1212    END DO
1213    !!     
1214  END SUBROUTINE Init_mask
1215  !
1216  !*****************************************************
1217  !   subroutine Init_Tmask(name,Grid)
1218  !*****************************************************
1219  !
1220  SUBROUTINE Init_Tmask(name,Grid,jpiglo,jpjglo)
1221    !
1222    USE io_netcdf
1223    !     
1224    CHARACTER(*) name
1225    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
1226    TYPE(Coordinates) :: Grid
1227    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
1228    !
1229    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
1230       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level)
1231    ELSE
1232       CALL read_ncdf_var('Bathy_level',name,Grid%Bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
1233    ENDIF
1234    !
1235    nx = SIZE(Grid%Bathy_level,1)
1236    ny = SIZE(Grid%Bathy_level,2) 
1237    !
1238    WRITE(*,*) 'Init masks in T points'   
1239    !     
1240    ALLOCATE( Grid%tmask(nx,ny,N) )
1241    !
1242    DO k = 1,N
1243       !     
1244       WHERE(Grid%Bathy_level(:,:) <= k-1 )   
1245          Grid%tmask(:,:,k) = 0.
1246       ELSEWHERE
1247          Grid%tmask(:,:,k) = 1.
1248       END WHERE
1249       !
1250    END DO
1251    !     
1252  END SUBROUTINE Init_Tmask
1253  !
1254  !*****************************************************
1255  !   subroutine Init_ssmask(name,Grid)
1256  !*****************************************************
1257  !
1258!  SUBROUTINE Init_ssmask(varname,filename,Grid,jpiglo,jpjglo)
1259!    !
1260!    USE io_netcdf
1261!    !     
1262!    CHARACTER(*) varname,filename
1263!    INTEGER :: nx,ny,k,i,j,jpiglo,jpjglo
1264!    TYPE(Coordinates) :: Grid
1265!    REAL*8, POINTER, DIMENSION(:,:) ::zwf => NULL()
1266!    !
1267!    IF(jpiglo == 1 .AND. jpjglo == 1) THEN
1268!       CALL read_ncdf_var(varname,filename,Grid%bathy_level)
1269!    ELSE
1270!       CALL read_ncdf_var(varname,filename,Grid%bathy_level,(/jpizoom,jpjzoom/),(/jpiglo,jpjglo/) )
1271!    ENDIF
1272!    !
1273!    nx = SIZE(Grid%bathy_level,1)
1274!    ny = SIZE(Grid%bathy_level,2)
1275!    !
1276!    WRITE(*,*) 'Init surface masks in T points'   
1277!    !     
1278!    ALLOCATE( Grid%ssmask(nx,ny), Grid%ssumask(nx,ny), Grid%ssvmask(nx,ny) )
1279!    !
1280!    WHERE(Grid%bathy_level(:,:) <= 0. )   
1281!       Grid%ssmask(:,:) = 0.
1282!    ELSEWHERE
1283!       Grid%ssmask(:,:) = 1.
1284!    END WHERE
1285!    !
1286!    Grid%ssumask(1:nx-1,:) = Grid%ssmask(1:nx-1,:)*Grid%ssmask(2:nx,:)
1287!    Grid%ssvmask(:,1:ny-1) = Grid%ssmask(:,1:ny-1)*Grid%ssmask(:,2:ny)
1288!    !
1289!    Grid%ssumask(nx,:) = Grid%ssmask(nx,:)
1290!    Grid%ssvmask(:,ny) = Grid%ssmask(:,ny)
1291!    !     
1292!  END SUBROUTINE Init_ssmask
1293  !
1294  !*****************************************************
1295  !   subroutine get_mask(name,Grid)
1296  !*****************************************************
1297  !
1298  SUBROUTINE get_mask(level,posvar,mask,filename)
1299    !
1300    USE io_netcdf
1301    !     
1302    CHARACTER(*) filename
1303    CHARACTER(*) posvar
1304    INTEGER :: level, nx, ny
1305    LOGICAL,DIMENSION(:,:),POINTER :: mask
1306    INTEGER,DIMENSION(:,:),POINTER :: maskT,maskU,maskV
1307    !     
1308    TYPE(Coordinates) :: Grid
1309    !
1310    CALL read_ncdf_var('Bathy_level',filename,Grid%Bathy_level)
1311    !
1312    nx = SIZE(Grid%Bathy_level,1)
1313    ny = SIZE(Grid%Bathy_level,2)
1314    ALLOCATE(maskT(nx,ny),mask(nx,ny))
1315    mask = .TRUE.
1316    !
1317    WHERE(Grid%Bathy_level(:,:) <= level-1 )   
1318       maskT(:,:) = 0
1319    ELSEWHERE
1320       maskT(:,:) = 1
1321    END WHERE
1322    !
1323    SELECT CASE(posvar)
1324       !
1325    CASE('T')
1326       !
1327       WHERE(maskT > 0)
1328          mask = .TRUE.
1329       ELSEWHERE
1330          mask = .FALSE.
1331       END WHERE
1332       DEALLOCATE(maskT)
1333       !
1334    CASE('U')
1335       !
1336       ALLOCATE(maskU(nx,ny))
1337       maskU(1:nx-1,:) = maskT(1:nx-1,:)*maskT(2:nx,:)
1338       maskU(nx,:) = maskT(nx,:)
1339       WHERE(maskU > 0)
1340          mask = .TRUE.
1341       ELSEWHERE
1342          mask = .FALSE.
1343       END WHERE
1344       DEALLOCATE(maskU,maskT)
1345       !
1346    CASE('V')
1347       !
1348       ALLOCATE(maskV(nx,ny))
1349       maskV(:,1:ny-1) = maskT(:,1:ny-1)*maskT(:,2:ny)
1350       maskV(:,ny) = maskT(:,ny)
1351       WHERE(maskT > 0)
1352          mask = .TRUE.
1353       ELSEWHERE
1354          mask = .FALSE.
1355       END WHERE
1356       DEALLOCATE(maskV,maskT)
1357       !
1358    END SELECT
1359    !     
1360  END SUBROUTINE get_mask
1361  !
1362  !
1363  !*****************************************************
1364  !   subroutine read_dimg_var(unit,irec,field)
1365  !*****************************************************
1366  !
1367  SUBROUTINE read_dimg_var(unit,irec,field,jpk)
1368    !     
1369    INTEGER :: unit,irec,jpk
1370    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
1371    INTEGER :: k
1372    !     
1373    DO k = 1,jpk
1374       READ(unit,REC=irec) field(:,:,k,1)
1375       irec = irec + 1
1376    ENDDO
1377    !
1378  END SUBROUTINE read_dimg_var
1379  !
1380  !
1381  !*****************************************************
1382  !   subroutine read_dimg_var(unit,irec,field)
1383  !*****************************************************
1384  !
1385  SUBROUTINE write_dimg_var(unit,irec,field,jpk)
1386    !     
1387    INTEGER :: unit,irec,jpk
1388    REAL*8,DIMENSION(:,:,:,:),POINTER :: field
1389    INTEGER :: k
1390    !     
1391    DO k = 1,jpk
1392       WRITE(unit,REC=irec) field(:,:,k,1)
1393       irec = irec + 1
1394    ENDDO
1395    !
1396  END SUBROUTINE write_dimg_var
1397
1398END MODULE agrif_readwrite
Note: See TracBrowser for help on using the repository browser.