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

Context Navigation

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

Last change on this file since 11206 was 11206, checked in by clem, 5 years ago
solve ticket #2280 by adding a new namelist parameter in the NESTING namelist to account for jperio of the parent grid
Property svn:keywords set to `Id`
File size: 53.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: