Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dombat.F90 @ 12101

Last change on this file since 12101 was 12101, checked in by mathiot, 4 years ago
merge ENHANCE-03_domcfg and ENHANCE-02_ISF_nemo
File size: 9.3 KB

Line
1	MODULE dombat
2
3	USE dom_oce ! ocean domain
4	! USE closea ! closed seas
5	!
6	USE in_out_manager ! I/O manager
7	USE iom ! I/O library
8	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
9	USE lib_mpp ! distributed memory computing library
10	USE agrif_modutil
11	USE bilinear_interp
12
13	IMPLICIT NONE
14	PRIVATE
15
16	PUBLIC dom_bat ! called by dom_zgr.F90
17
18	CONTAINS
19
20	SUBROUTINE dom_bat
21
22	INTEGER :: inum, isize, jsize, id, ji, jj
23	INTEGER :: tabdim1, tabdim2, nxhr, nyhr, nxyhr
24	INTEGER, DIMENSION(2) :: ddims
25	INTEGER, DIMENSION(3) :: status
26	INTEGER, DIMENSION(:,:), ALLOCATABLE :: trouble_points , vardep,mask_oce
27	INTEGER :: iimin,iimax,jjmin,jjmax
28	INTEGER, DIMENSION(1) :: i_min,i_max
29	INTEGER, DIMENSION(1) ::j_min,j_max
30	REAL(wp), DIMENSION(jpi,jpj) :: myglamf
31	INTEGER,DIMENSION(:) ,POINTER :: src_add,dst_add => NULL()
32	REAL(wp), DIMENSION(:) ,ALLOCATABLE :: vardep1d, lon_new1D,lat_new1D
33	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bathy_new, lat_new, lon_new, bathy_test
34	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: coarselon, coarselat, coarsebathy
35	REAL(wp) :: Cell_lonmin, Cell_lonmax, Cell_latmin, Cell_latmax
36	LOGICAL :: identical_grids
37	LOGICAL, DIMENSION(:,:), ALLOCATABLE :: masksrc
38	REAL*8, DIMENSION(:,:),POINTER :: matrix,interpdata => NULL()
39	LOGICAL :: lonlat_2D
40
41	CHARACTER(32) :: bathyfile, bathyname, lonname,latname
42
43	lonlat_2D=.false.
44
45	bathyfile=TRIM(cn_topo)
46	bathyname=TRIM(cn_bath)
47	lonname=TRIM(cn_lon)
48	latname=TRIM(cn_lat)
49
50	CALL iom_open( bathyfile, inum, lagrif=.FALSE. )
51	id = iom_varid( inum, bathyname, ddims )
52
53	status=-1
54	ALLOCATE(lon_new (ddims(1),ddims(2)), STAT=status(1))
55	ALLOCATE(lat_new (ddims(1),ddims(2)), STAT=status(2))
56	ALLOCATE(bathy_new(ddims(1),ddims(2)), STAT=status(3))
57	IF( sum(status) /= 0 ) CALL ctl_stop( 'STOP', 'dom_bat : unable to allocate arrays' )
58
59	IF (lonlat_2D) THEN
60	CALL iom_get ( inum, jpdom_unknown, lonname, lon_new )
61	CALL iom_get ( inum, jpdom_unknown, latname, lat_new )
62	ELSE
63	ALLOCATE(lon_new1D(ddims(1)), lat_new1D(ddims(2)))
64	CALL iom_get ( inum, jpdom_unknown, lonname, lon_new1D )
65	CALL iom_get ( inum, jpdom_unknown, latname, lat_new1D )
66	DO ji=1, ddims(1)
67	lon_new(ji,:)=lon_new1D(ji)
68	ENDDO
69	DO ji=1, ddims(2)
70	lat_new(:,ji)=lat_new1D(ji)
71	ENDDO
72	ENDIF
73	CALL iom_get ( inum, jpdom_unknown, bathyname, bathy_new )
74	CALL iom_close (inum)
75	WHERE (bathy_new > 0.)
76	bathy_new=0.
77	ENDWHERE
78	bathy_new=-bathy_new
79
80	! Eventually add here a pre-selection of the area (coarselon/lat)
81	i_min=10 ; j_min=10
82	i_max= ddims(1)-10 ; j_max=ddims(2)-10
83
84	tabdim1 = i_max(1) - i_min(1) + 1
85	tabdim2 = j_max(1) - j_min(1) + 1
86	!
87
88	ALLOCATE(coarselon(tabdim1,tabdim2))
89	ALLOCATE(coarselat(tabdim1,tabdim2))
90	ALLOCATE(coarsebathy(tabdim1,tabdim2))
91	!
92	WHERE( lon_new < 0. ) lon_new = lon_new + 360.
93	myglamf=glamf
94	WHERE( myglamf < 0. ) myglamf = myglamf + 360.
95
96	coarselat(:,:) = lat_new (i_min(1):i_max(1), j_min(1):j_max(1))
97	coarselon (:,:) = lon_new (i_min(1):i_max(1), j_min(1):j_max(1))
98	coarsebathy(:,:) = bathy_new(i_min(1):i_max(1), j_min(1):j_max(1))
99
100	IF( nn_interp == 0 .OR. nn_interp == 1 ) THEN ! arithmetic or median averages
101	! ! -----------------------------
102	ALLOCATE(trouble_points(jpi,jpj))
103	trouble_points(:,:) = 0
104	!
105	DO jj = 2, jpj
106	DO ji = 2, jpi
107	!
108	! fine grid cell extension
109	Cell_lonmin = MIN(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj))
110	Cell_lonmax = MAX(myglamf(ji-1,jj-1),myglamf(ji,jj-1),myglamf(ji,jj),myglamf(ji-1,jj))
111	Cell_latmin = MIN(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj))
112	Cell_latmax = MAX(gphif(ji-1,jj-1),gphif(ji,jj-1),gphif(ji,jj),gphif(ji-1,jj))
113	!
114	! look for points in G0 (bathy dataset) contained in the fine grid cells
115	iimin = 1
116	DO WHILE( coarselon(iimin,1) < Cell_lonmin )
117	iimin = iimin + 1
118	ENDDO
119	!
120	jjmin = 1
121	DO WHILE( coarselat(iimin,jjmin) < Cell_latmin )
122	jjmin = jjmin + 1
123	ENDDO
124	!
125	iimax = iimin
126	DO WHILE( coarselon(iimax,1) <= Cell_lonmax )
127	iimax = iimax + 1
128	iimax = MIN( iimax,SIZE(coarsebathy,1))
129	ENDDO
130	!
131	jjmax = jjmin
132	DO WHILE( coarselat(iimax,jjmax) <= Cell_latmax )
133	jjmax = jjmax + 1
134	jjmax = MIN( jjmax,SIZE(coarsebathy,2))
135	ENDDO
136	!
137	IF( .NOT. Agrif_Root() ) THEN
138	iimax = iimax-1
139	jjmax = jjmax-1
140	ELSE
141	iimax = MAX(iimin,iimax-1)
142	jjmax = MAX(jjmin,jjmax-1)
143	ENDIF
144	!
145	iimin = MAX( iimin,1 )
146	jjmin = MAX( jjmin,1 )
147	iimax = MIN( iimax,SIZE(coarsebathy,1))
148	jjmax = MIN( jjmax,SIZE(coarsebathy,2))
149
150	nxhr = iimax - iimin + 1
151	nyhr = jjmax - jjmin + 1
152
153
154	IF( nxhr == 0 .OR. nyhr == 0 ) THEN
155	!
156	trouble_points(ji,jj) = 1
157	!
158	ELSE
159	!
160	ALLOCATE( vardep(nxhr,nyhr), mask_oce(nxhr,nyhr) )
161	vardep(:,:) = coarsebathy(iimin:iimax,jjmin:jjmax)
162	!
163	WHERE( vardep(:,:) .GT. 0. ) ; mask_oce = 1 ;
164	ELSEWHERE ; mask_oce = 0 ;
165	ENDWHERE
166	!
167	nxyhr = nxhr*nyhr
168	IF( SUM(mask_oce) < 0.5*(nxyhr) ) THEN ! if more than half of the points are on land then bathy fine = 0
169	bathy(ji,jj) = 0.
170	ELSE
171	IF( nn_interp == 0 ) THEN ! Arithmetic average
172	bathy(ji,jj) = SUM( vardep(:,:) * mask_oce(:,:) ) / SUM( mask_oce(:,:) )
173	ELSE ! Median average
174	ALLOCATE(vardep1d(nxyhr))
175	vardep1d = RESHAPE(vardep,(/ nxyhr /) )
176	!!CALL ssort(vardep1d,nxyhr)
177	CALL quicksort(vardep1d,1,nxyhr)
178	!
179	! Calculate median
180	IF (MOD(nxyhr,2) .NE. 0) THEN
181	bathy(ji,jj) = vardep1d( nxyhr/2 + 1 )
182	ELSE
183	bathy(ji,jj) = 0.5 * ( vardep1d(nxyhr/2) + vardep1d(nxyhr/2+1) )
184	END IF
185	DEALLOCATE(vardep1d)
186	ENDIF
187	ENDIF
188	DEALLOCATE (mask_oce,vardep)
189	!
190	ENDIF
191	ENDDO
192	ENDDO
193
194	IF( SUM( trouble_points ) > 0 ) THEN
195	PRINT*,'too much empty cells, proceed to bilinear interpolation'
196	nn_interp = 2
197	stop
198	ENDIF
199	ENDIF
200
201	#undef MYTEST
202	#ifdef MYTEST
203	IF( nn_interp == 2) THEN ! Bilinear interpolation
204	! ! -----------------------------
205	identical_grids = .FALSE.
206
207	IF( SIZE(coarselat,1) == jpi .AND. SIZE(coarselat,2) == jpj .AND. &
208	SIZE(coarselon,1) == jpj .AND. SIZE(coarselon,2) == jpj ) THEN
209	IF( MAXVAL( ABS(coarselat(:,:)- gphit(:,:)) ) < 0.0001 .AND. &
210	MAXVAL( ABS(coarselon(:,:)- glamt(:,:)) ) < 0.0001 ) THEN
211	PRINT*,'same grid between ', cn_topo,' and child domain'
212	bathy = bathy_new
213	identical_grids = .TRUE.
214	ENDIF
215	ENDIF
216
217	IF( .NOT. identical_grids ) THEN
218
219	ALLOCATE(masksrc(SIZE(bathy_new,1),SIZE(bathy_new,2)))
220	ALLOCATE(bathy_test(jpi,jpj))
221	!
222	!Where(G0%bathy_meter.le.0.00001)
223	! masksrc = .false.
224	!ElseWhere
225	masksrc = .TRUE.
226	!End where
227	!
228	! compute remapping matrix thanks to SCRIP package
229	CALL get_remap_matrix(coarselat,gphit,coarselon,glamt,masksrc,matrix,src_add,dst_add)
230	CALL make_remap(bathy_new,bathy_test,jpi,jpj,matrix,src_add,dst_add)
231	!
232	bathy = bathy_test
233	!
234	DEALLOCATE(masksrc)
235	DEALLOCATE(bathy_test)
236
237	ENDIF
238	!
239	ENDIF
240	#endif
241	END SUBROUTINE dom_bat
242
243	END MODULE dombat

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

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

Context Navigation

source: utils/tools_UKMO_MERGE_2019/DOMAINcfg/src/dombat.F90 @ 12101

Download in other formats: