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

Last change on this file since 10383 was 10383, checked in by clem, 5 years ago

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