Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

agrif_interpolation.f90 @ 14974

Last change on this file since 14974 was 14974, checked in by clem, 3 years ago
NESTING TOOLS 4.0-HEAD: solve ticket #2687
Property svn:keywords set to `Id`
File size: 48.5 KB

Line
1	!************************************************************************
2	! Fortran 95 OPA Nesting tools *
3	! *
4	! Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr) *
5	! *
6	!************************************************************************
7	!
8	MODULE agrif_interpolation
9	!
10	USE agrif_types
11	USE io_netcdf
12	USE bicubic_interp
13	USE bilinear_interp
14	USE agrif_extrapolation
15	USE agrif_readwrite
16	!
17	!
18	!************************************************************************
19	! *
20	! MODULE AGRIF_INTERPOLATION *
21	! *
22	! module containing subroutine used for : *
23	! - Forcing data interpolation *
24	! - Parent to Child coordinates interpolation *
25	! *
26	!************************************************************************
27	!
28	!
29	CONTAINS
30	!
31	!
32	!****************************************************************
33	! subroutine agrif_interp *
34	! *
35	! subroutine to interpolate coordinates *
36	! *
37	! - input : *
38	! tabin : coarse grid coordinate variable *
39	! typevar : position of interpolated variable on cells *
40	! *
41	! - ouput : *
42	! tabout : coordinate variable interpolated on fine grid *
43	! *
44	!****************************************************************
45	!
46	SUBROUTINE agrif_interp(tabin,tabout,typevar)
47	!
48	REAL*8,DIMENSION(:,:) :: tabin,tabout
49	CHARACTER(*) :: typevar
50	REAL*8 :: cff1
51	INTEGER :: ii,jj
52	REAL*8,DIMENSION(:,:),ALLOCATABLE :: tabouttemp
53	!
54	!
55	!
56	IF( ln_agrif_domain ) THEN
57	CALL agrif_base_interp2(tabin,tabout,imin,jmin,typevar)
58	ELSE
59	CALL agrif_base_interp3(tabin,tabout,imin,jmin,typevar)
60	ENDIF
61
62	!
63	END SUBROUTINE agrif_interp
64	!
65	!
66	SUBROUTINE agrif_base_interp2(tabin,tabout,i_min,j_min,typevar)
67	!
68	IMPLICIT NONE
69	REAL*8,DIMENSION(:,:) :: tabin,tabout
70	REAL*8 :: cff1
71	INTEGER :: i_min,j_min
72	INTEGER :: ii,jj,i,j
73	CHARACTER(*) :: typevar
74	REAL*8,DIMENSION(:),ALLOCATABLE :: xc,yc,xf,yf
75	REAL*8,DIMENSION(:,:),ALLOCATABLE :: tabtemp
76	REAL*8 :: dxc,dyc,dxf,dyf
77	REAL*8 :: decalxc,decalyc,decalxf,decalyf
78
79	INTEGER :: ptx,pty
80	REAL*8 :: val(4),xval(4)
81
82	INTEGER :: nxc,nyc,nxf,nyf
83	REAL :: xmin,ymin,x,y
84	INTEGER :: ic,jc,itemp,jtemp
85
86	nxc = SIZE(tabin,1)
87	nyc = SIZE(tabin,2)
88
89	nxf = SIZE(tabout,1)
90	nyf = SIZE(tabout,2)
91	!
92
93	ALLOCATE(xc(nxc))
94	ALLOCATE(yc(nyc))
95	ALLOCATE(xf(nxf))
96	ALLOCATE(yf(nyf))
97
98	ALLOCATE(tabtemp(nxc,nyf))
99
100	dxc = 1.
101	dyc = 1.
102	dxf = 1./REAL(irafx)
103	dyf = 1./REAL(irafy)
104
105	IF (typevar .EQ. 'F') THEN
106	ptx = 1 + nbghostcellsfine
107	pty = 1 + nbghostcellsfine
108	decalxc = 0.
109	decalyc = 0.
110	decalxf = 0.
111	decalyf = 0.
112	ELSEIF (typevar .EQ. 'T') THEN
113	ptx = 2 + nbghostcellsfine
114	pty = 2 + nbghostcellsfine
115	decalxc = dxc/2.
116	decalyc = dyc/2.
117	decalxf = dxf/2.
118	decalyf = dyf/2.
119	ELSEIF (typevar .EQ. 'U') THEN
120	ptx = 1 + nbghostcellsfine
121	pty = 2 + nbghostcellsfine
122	decalxc = 0.
123	decalyc = dyc/2.
124	decalxf = 0.
125	decalyf = dyf/2.
126	ELSEIF (typevar .EQ. 'V') THEN
127	ptx = 2 + nbghostcellsfine
128	pty = 1 + nbghostcellsfine
129	decalxc = dxc/2.
130	decalyc = 0.
131	decalxf = dxf/2.
132	decalyf = 0.
133	ENDIF
134
135	DO i=1,nxc
136	xc(i) = 0. + (i-ptx) * dxc + decalxc
137	ENDDO
138
139	DO j=1,nyc
140	yc(j) = 0. + (j-pty) * dyc + decalyc
141	ENDDO
142
143
144	xmin = (i_min - 1) * dxc
145	ymin = (j_min - 1) * dyc
146
147	DO i=1,nxf
148	xf(i) = xmin + (i-ptx) * dxf + decalxf
149	ENDDO
150
151	DO j=1,nyf
152	yf(j) = ymin + (j-pty) * dyf + decalyf
153	ENDDO
154
155
156	DO j = 1,nyf
157	DO i = 1,nxc
158	x = xc(i)
159	y = yf(j)
160
161	ic = ptx + NINT((x-0.-decalxc)/dxc)
162	jc = pty + agrif_int((y-0.-decalyc)/dyc)
163
164	jc = jc - 1
165
166	CALL polint(yc(jc:jc+3),tabin(ic,jc:jc+3),4,y,tabtemp(i,j))
167	ENDDO
168	ENDDO
169
170	DO j = 1,nyf
171	DO i = 1,nxf
172	x = xf(i)
173	y = yf(j)
174
175	itemp = ptx + agrif_int((x-0.-decalxc)/dxc)
176	jtemp = pty + NINT((y-ymin-decalyf)/dyf)
177
178	itemp = itemp - 1
179
180	val = tabtemp(itemp:itemp+3,jtemp)
181	xval = xc(itemp:itemp+3)
182	CALL polint(xval,val,4,x,tabout(i,j))
183
184	ENDDO
185	ENDDO
186
187	DEALLOCATE(xc,yc,xf,yf,tabtemp)
188
189
190	END SUBROUTINE agrif_base_interp2
191	!
192	!
193	SUBROUTINE agrif_base_interp3(tabin,tabout,i_min,j_min,typevar)
194	!
195	IMPLICIT NONE
196	REAL*8,DIMENSION(:,:) :: tabin,tabout
197	INTEGER :: i_min,j_min
198	CHARACTER(*) :: typevar
199
200	INTEGER :: nxf,nyf,nxc,nyc,zx,zy
201	INTEGER :: ji,jj,jif,jjf,jic,jjc,jic1,jjc1,jdecx,jdecy
202	REAL*8 :: Ax, Bx, Ay, By
203
204	nxf = SIZE(tabout,1)
205	nyf = SIZE(tabout,2)
206
207	nxc = SIZE(tabin,1)
208	nyc = SIZE(tabin,2)
209
210	SELECT CASE(typevar)
211	CASE('T')
212	IF(MOD(irafx,2)==1) THEN ! odd
213	zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
214	ELSE ! even
215	zx = 2 ; zy = 2 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
216	ENDIF
217	CASE('U')
218	IF(MOD(irafx,2)==1) THEN ! odd
219	zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
220	ELSE ! even
221	zx = 1 ; zy = 2 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
222	ENDIF
223	CASE('V')
224	IF(MOD(irafx,2)==1) THEN ! odd
225	zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
226	ELSE ! even
227	zx = 2 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
228	ENDIF
229	CASE('F')
230	IF(MOD(irafx,2)==1) THEN ! odd
231	zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
232	ELSE ! even
233	zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
234	ENDIF
235	END SELECT
236
237
238	DO jj = 1, nyf
239
240	jjf = jj - jdecy
241	jjc = j_min + FLOOR((jjf-1.) / irafy)
242
243	DO ji = 1, nxf
244
245	jif = ji - jdecx
246	jic = i_min + FLOOR((jif-1.) / irafx)
247
248	Bx = MOD( zxjif-1, zxirafx ) / REAL(zx*irafx)
249	By = MOD( zyjjf-1, zyirafy ) / REAL(zy*irafy)
250	Ax = 1. - Bx
251	Ay = 1. - By
252
253	jic1 = MIN( nxc, jic+1 ) ! avoid out of bounds for tabin below
254	jjc1 = MIN( nyc, jjc+1 ) ! --
255
256	tabout(ji,jj) = ( Bx * tabin(jic1,jjc ) + Ax * tabin(jic,jjc ) ) * Ay + &
257	& ( Bx * tabin(jic1,jjc1) + Ax * tabin(jic,jjc1) ) * By
258	END DO
259	END DO
260
261	!
262	END SUBROUTINE agrif_base_interp3
263
264	!
265	SUBROUTINE polint(xin,valin,n,x,val)
266	IMPLICIT NONE
267	INTEGER n
268	REAL*8 xin(n), valin(n)
269	REAL*8 x,val
270
271	INTEGER ns,i,m
272	REAL *8 dif,dift
273	REAL*8 c(n),d(n),ho,hp,w,den,dy
274
275	ns = 1
276	dif = ABS(x-xin(1))
277
278	DO i=1,n
279	dift = ABS(x-xin(i))
280	IF (dift < dif) THEN
281	ns = i
282	dif = dift
283	ENDIF
284	c(i) = valin(i)
285	d(i) = valin(i)
286	ENDDO
287
288	val = valin(ns)
289	ns = ns - 1
290
291	DO m=1,n-1
292	DO i=1,n-m
293	ho = xin(i)-x
294	hp = xin(i+m)-x
295	w=c(i+1)-d(i)
296	den = w/(ho-hp)
297	d(i) = hp * den
298	c(i) = ho * den
299	ENDDO
300	IF (2*ns < (n-m)) THEN
301	dy = c(ns+1)
302	ELSE
303	dy = d(ns)
304	ns = ns - 1
305	ENDIF
306	val = val + dy
307	ENDDO
308	END SUBROUTINE polint
309	!
310	!
311	SUBROUTINE polcoe(xin,valin,n,cof)
312	IMPLICIT NONE
313	INTEGER n
314	REAL*8 xin(n),valin(n),cof(n)
315
316
317	INTEGER i,j,k
318	REAL*8 b,ff,phi,s(n)
319
320	s = 0.
321	cof = 0.
322
323	s(n)=-xin(1)
324	DO i=2,n
325	DO j=n+1-i,n-1
326	s(j) =s(j) -xin(i)*s(j+1)
327	ENDDO
328	s(n)=s(n)-xin(i)
329	ENDDO
330
331	DO j=1,n
332	phi=n
333	DO k=n-1,1,-1
334	phi = ks(k+1)+xin(j)phi
335	ENDDO
336	ff=valin(j)/phi
337	b=1.
338	DO k=n,1,-1
339	cof(k)=cof(k)+b*ff
340	b=s(k)+xin(j)*b
341	ENDDO
342	ENDDO
343
344	RETURN
345	END SUBROUTINE polcoe
346	!
347
348	!****************************************************************
349	! subroutine Correctforconservation *
350	! *
351	! Conservation on domain boundaries ( over 2 coarse grid cells) *
352	! *
353	!****************************************************************
354	!
355	!
356	SUBROUTINE Correctforconservation(tabcoarse,tabfine,e1parent,e2parent,e1,e2,nxfin,nyfin,posvar,i_min,j_min)
357	IMPLICIT NONE
358	INTEGER :: nxfin,nyfin
359	REAL*8,DIMENSION(:,:) :: tabcoarse,tabfine,e1parent,e2parent,e1,e2
360	CHARACTER*1 :: posvar
361	INTEGER :: i_min,j_min,diff
362	INTEGER ji,jj,ipt,jpt,i,j
363	INTEGER ind1,ind2,ind3,ind4,ind5,ind6
364	REAL*8 cff1,cff2,cff3
365	!
366	diff = 0
367	IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
368	!
369	ind1 = 2 + CEILING(irafx/2.0) + diff
370	ind2 = nxfin-(2-1)-CEILING(irafx/2.0)
371	ind3 = 2 + CEILING(irafy/2.0) + diff
372	ind4 = nyfin-(2-1)-CEILING(irafy/2.0)
373	ind5 = nxfin - 1 - irafx - CEILING(irafx/2.0)
374	ind6 = nyfin - 1 - irafy - CEILING(irafy/2.0)
375	!
376	IF (posvar.EQ.'T') THEN
377	!
378	PRINT *,'correction'
379	!
380	DO ji=ind1,ind1+irafx,irafx
381	!
382	ipt=i_min+(3-1)+(ji-ind1)/irafx
383	!
384	DO jj=ind3,ind4,irafy
385	!
386	jpt=j_min+(3-1)+(jj-ind3)/irafy
387
388	cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
389	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
390
391	cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
392	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
393	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
394
395	cff3=e1parent(ipt,jpt)e2parent(ipt,jpt)tabcoarse(ipt,jpt)
396
397	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)= &
398	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1
399
400	END DO
401
402	ENDDO
403	!****
404	DO ji=ind5,ind5+irafx,irafx
405	!
406	ipt=i_min+(3-1)+(ji-ind1)/irafx
407
408	DO jj=ind3,ind4,irafy
409	jpt=j_min+(3-1)+(jj-ind3)/irafy
410
411	cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
412	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
413
414	cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
415	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
416	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
417
418	cff3=e1parent(ipt,jpt)e2parent(ipt,jpt)tabcoarse(ipt,jpt)
419
420	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)= &
421	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1
422
423	END DO
424	ENDDO
425
426	DO jj=ind3,ind3+irafy,irafy
427	!
428	jpt=j_min+(3-1)+(jj-ind3)/irafy
429	!
430	DO ji=ind1,ind2,irafx
431	!
432	ipt=i_min+(3-1)+(ji-ind1)/irafx
433
434	cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
435	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
436
437	cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
438	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
439	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
440
441	cff3=e1parent(ipt,jpt)e2parent(ipt,jpt)tabcoarse(ipt,jpt)
442
443	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)= &
444	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1
445
446	END DO
447	ENDDO
448	!****
449	DO jj=ind6,ind6+irafy,irafy
450	!
451	jpt=j_min+(3-1)+(jj-ind3)/irafy
452	!
453	DO ji=ind1,ind2,irafx
454	ipt=i_min+(3-1)+(ji-ind1)/irafx
455
456	cff1=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
457	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
458
459	cff2=SUM(e1(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
460	e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)* &
461	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
462
463	cff3=e1parent(ipt,jpt)e2parent(ipt,jpt)tabcoarse(ipt,jpt)
464
465	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)= &
466	tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)+(cff3-cff2)/cff1
467
468	END DO
469	ENDDO
470	!
471	!
472	ENDIF
473	RETURN
474	END SUBROUTINE Correctforconservation
475
476	!
477	!****************************************************************
478	! subroutine Interp_Extrap_var *
479	! *
480	! Interpolation and Extrapolation for temperature and salinity *
481	! *
482	! -input : *
483	! filename : file containing variable to interpolate *
484	! *
485	!****************************************************************
486	!
487	SUBROUTINE Interp_Extrap_var(filename, cd_type)
488	!
489	IMPLICIT NONE
490	!
491	! variables declaration
492	!
493	CHARACTER(len=80),INTENT(in) :: filename
494	CHARACTER(len=1 ), OPTIONAL, INTENT(in) :: cd_type
495
496	CHARACTER*100 :: interp_type
497	CHARACTER*100 :: Child_file,Childcoordinates
498	CHARACTER*100 :: varname,childbathy
499	CHARACTER*1 :: posvar
500	CHARACTER*20,DIMENSION(:),POINTER :: Ncdf_varname
501	CHARACTER(len=20),DIMENSION(4) :: dimnames
502	!
503	REAL*8, POINTER, DIMENSION(:,:) :: lonParent,latParent => NULL()
504	REAL*8, POINTER, DIMENSION(:,:) :: lonChild,latChild,latlon_temp => NULL()
505	REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL()
506	REAL*8, POINTER, DIMENSION(:,:,:) :: tabinterp3d,tabvar3d,tabvar01,tabvar02,tabvar03 => NULL()
507	REAL*8, POINTER, DIMENSION(:) :: timedepth_temp,depth => NULL()
508	REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL()
509	INTEGER,DIMENSION(:),POINTER :: src_add,dst_add => NULL()
510	INTEGER, POINTER, DIMENSION(:) :: tabtemp1DInt => NULL()
511	REAL*8, POINTER, DIMENSION(:) :: nav_lev => NULL()
512	!
513	LOGICAL, DIMENSION(:,:,:), POINTER :: detected_pts => NULL()
514	REAL*8,DIMENSION(:,:,:),POINTER :: mask => NULL()
515	LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
516	LOGICAL :: Interpolation,conservation,Pacifique,Extrapolation,land_level
517	!
518	INTEGER :: deptht,time,i,status,ncid,t,ii,j,nb,numlon,numlat,unlimdimid
519
520	!
521	TYPE(Coordinates) :: G0,G1
522	!
523	status = nf90_open(TRIM(filename),NF90_NOWRITE,ncid)
524	IF (status/=nf90_noerr) THEN
525	WRITE(,)"unable to open netcdf file : ",TRIM(filename)
526	STOP
527	ENDIF
528	!
529	!*****************
530	!If coarse grid is masked possibility to activate an extrapolation process
531	!
532	Extrapolation = .FALSE.
533	PRINT*,'DEBUT INTERP_EXTRAP_VAR'
534	!
535	!*****************
536	!
537	! check grid position
538	!
539	IF( PRESENT(cd_type) ) THEN
540	posvar = cd_type
541	ELSE
542	posvar = 'T'
543	ENDIF
544
545	!
546	! read dimensions in netcdf file
547	!
548	CALL Read_Ncdf_dim('x',filename,numlon)
549	CALL Read_Ncdf_dim('y',filename,numlat)
550	IF ( Dims_Existence( 'deptht' , filename ) ) THEN
551	CALL Read_Ncdf_dim('deptht',filename,deptht)
552	ELSE IF ( Dims_Existence( 'depthu' , filename ) ) THEN
553	CALL Read_Ncdf_dim('depthu',filename,deptht)
554	ELSE IF ( Dims_Existence( 'depthv' , filename ) ) THEN
555	CALL Read_Ncdf_dim('depthv',filename,deptht)
556	ELSE IF ( Dims_Existence( 'depthw' , filename ) ) THEN
557	CALL Read_Ncdf_dim('depthw',filename,deptht)
558	ELSE IF ( Dims_Existence( 'z' , filename ) ) THEN
559	CALL Read_Ncdf_dim('z',filename,deptht)
560	ELSE
561	deptht = N
562	ENDIF
563	IF ( Dims_Existence( 'time_counter' , filename ) ) THEN
564	CALL Read_Ncdf_dim('time_counter',filename,time)
565	! check that time_counter is the unlimited dim
566	status = nf90_inquire(ncid, unlimiteddimid = unlimdimid)
567	IF ( unlimdimid == -1 ) THEN
568	WRITE(,)"time_counter should be the unlimited dimension in this file : ",filename
569	WRITE(,)" use nco command: ncks -O --mk_rec_dmn time_counter ",filename," ",filename
570	STOP
571	ENDIF
572	ELSEIF( Dims_Existence( 'time' , filename ) ) THEN
573	CALL Read_Ncdf_dim('time',filename,time)
574	! check that time is the unlimited dim
575	status = nf90_inquire(ncid, unlimiteddimid = unlimdimid)
576	IF ( unlimdimid == -1 ) THEN
577	WRITE(,)"time should be the unlimited dimension in this file : ",filename
578	WRITE(,)" use nco command: ncks -O --mk_rec_dmn time ",filename," ",filename
579	STOP
580	ENDIF
581	ELSE
582	time = 0
583	ENDIF
584
585	!
586	! retrieve netcdf variable name
587	!
588	CALL Read_Ncdf_VarName(filename,Ncdf_varname)
589	!
590	! define name of child grid file
591	!
592	CALL set_child_name(filename,Child_file)
593	CALL set_child_name(parent_coordinate_file,Childcoordinates)
594	CALL set_child_name(parent_bathy_level,childbathy)
595	WRITE(,) 'Child grid file name = ',TRIM(Child_file)
596	!
597	! create this file
598	!
599	status = nf90_create(Child_file,NF90_WRITE,ncid)
600	status = nf90_close(ncid)
601	!
602	! read coordinates of both domains
603	!
604	status = Read_Coordinates(parent_coordinate_file,G0)
605	status = Read_Coordinates(Childcoordinates,G1,Pacifique)
606	!
607	! check consistency of informations read in namelist
608	!
609	IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
610	imax <= imin .OR. jmax <= jmin ) THEN
611	WRITE(,) 'ERROR ***** bad child grid definition ...'
612	WRITE(,) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'
613	WRITE(,) ' '
614	STOP
615	ENDIF
616	!
617	IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
618	WRITE(,) 'ERROR ***** bad child coordinates or bathy file ...'
619	WRITE(,) ' '
620	WRITE(,) 'please check that child coordinates file and child bathymetry file'
621	WRITE(,) 'has been created with the current namelist '
622	WRITE(,) ' '
623	STOP
624	ENDIF
625	!
626
627	!
628	! Initialization of T-mask thanks to bathymetry
629	!
630	IF( Extrapolation ) THEN
631
632	WRITE(,) 'mask initialisation on coarse and fine grids'
633	!
634	ALLOCATE(mask(numlon,numlat,N))
635	CALL Init_mask(childbathy,G1,1,1)
636	CALL Init_mask(parent_bathy_level,G0,1,1)
637	!
638	ENDIF
639
640	!
641	! select coordinates to use according to variable position
642	!
643	ALLOCATE(lonParent(numlon,numlat),latParent(numlon,numlat))
644	ALLOCATE(lonChild(nxfin,nyfin),latChild(nxfin,nyfin))
645
646	SELECT CASE(posvar)
647	CASE('T')
648	lonParent = G0%glamt
649	latParent = G0%gphit
650	lonChild = G1%glamt
651	latChild = G1%gphit
652	mask = G1%tmask
653	CASE('U')
654	lonParent = G0%glamu
655	latParent = G0%gphiu
656	lonChild = G1%glamu
657	latChild = G1%gphiu
658	mask = G1%umask
659	CASE('V')
660	lonParent = G0%glamv
661	latParent = G0%gphiv
662	lonChild = G1%glamv
663	latChild = G1%gphiv
664	mask = G1%vmask
665	END SELECT
666
667	DEALLOCATE(G0%glamu,G0%glamv,G0%gphiu,G0%gphiv)
668	DEALLOCATE(G1%glamu,G1%glamv,G1%gphiu,G1%gphiv)
669	DEALLOCATE(G1%glamt,G1%gphit,G0%glamt,G0%gphit)
670
671	!
672	!longitude modification if child domain covers Pacific ocean area
673	!
674	IF( lonChild(1,1) > lonChild(nxfin,nyfin) ) THEN
675	Pacifique = .TRUE.
676	WHERE( lonChild < 0 )
677	lonChild = lonChild + 360.
678	END WHERE
679	WHERE( lonParent < 0 )
680	lonParent = lonParent + 360.
681	END WHERE
682	ENDIF
683
684	!
685	!
686	! dimensions initialization
687	!
688	CALL Write_Ncdf_dim('x',Child_file,nxfin)
689	CALL Write_Ncdf_dim('y',Child_file,nyfin)
690	IF ( Dims_Existence( 'deptht' , filename ) ) CALL Write_Ncdf_dim('deptht',Child_file,deptht)
691	IF ( Dims_Existence( 'depthu' , filename ) ) CALL Write_Ncdf_dim('depthu',Child_file,deptht)
692	IF ( Dims_Existence( 'depthv' , filename ) ) CALL Write_Ncdf_dim('depthv',Child_file,deptht)
693	IF ( Dims_Existence( 'depthw' , filename ) ) CALL Write_Ncdf_dim('depthw',Child_file,deptht)
694	IF ( Dims_Existence( 'z' , filename ) ) CALL Write_Ncdf_dim('z' ,Child_file,deptht)
695	IF ( Dims_Existence( 'time_counter' , filename ) ) CALL Write_Ncdf_dim('time_counter',Child_file,time)
696	IF ( Dims_Existence( 'time' , filename ) ) CALL Write_Ncdf_dim('time' ,Child_file,time)
697
698	IF( deptht .NE. 1 .AND. deptht .NE. N ) THEN
699	WRITE(,) '***'
700	WRITE(,) 'Number of vertical levels doesn t match between namelist'
701	WRITE(,) 'and forcing file ',TRIM(filename)
702	WRITE(,) 'Please check the values in namelist file'
703	WRITE(,) 'N = ',N
704	WRITE(,) 'deptht = ',deptht
705	STOP
706	ENDIF
707	!
708	!
709	DO i = 1,SIZE(Ncdf_varname)
710	!
711	!
712	! ****************************LOOP ON VARIABLE NAMES*****************************************
713	!
714	!
715	SELECT CASE (TRIM(Ncdf_varname(i)))
716	!
717	!copy nav_lon from child coordinates to output file
718	CASE('nav_lon')
719	ALLOCATE(latlon_temp(nxfin,nyfin))
720	CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),latlon_temp)
721	CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,latlon_temp,'float')
722	CALL Copy_Ncdf_att('nav_lon',TRIM(filename),Child_file, &
723	MINVAL(latlon_temp),MAXVAL(latlon_temp))
724	DEALLOCATE(latlon_temp)
725	varname = TRIM(Ncdf_varname(i))
726	Interpolation = .FALSE.
727	!
728	!copy nav_lat from child coordinates to output file
729	CASE('nav_lat')
730	ALLOCATE(latlon_temp(nxfin,nyfin))
731	CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),latlon_temp)
732	CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,latlon_temp,'float')
733	CALL Copy_Ncdf_att('nav_lat',TRIM(filename),Child_file, &
734	MINVAL(latlon_temp),MAXVAL(latlon_temp))
735	DEALLOCATE(latlon_temp)
736	varname = TRIM(Ncdf_varname(i))
737	Interpolation = .FALSE.
738	!
739	!copy nav_lev from restart_file to output file
740	!
741	CASE('nav_lev')
742
743	WRITE(,) 'copy nav_lev'
744	CALL Read_Ncdf_var('nav_lev',filename,nav_lev)
745	IF(.NOT. dimg ) THEN
746	CALL Write_Ncdf_var('nav_lev','z',Child_file,nav_lev,'float')
747	CALL Copy_Ncdf_att('nav_lev',filename,Child_file)
748	ENDIF
749	DEALLOCATE(nav_lev)
750	Interpolation = .FALSE.
751	!
752	!copy time_counter from input file to output file
753	CASE('time_counter')
754	ALLOCATE(timedepth_temp(time))
755	CALL Read_Ncdf_var('time_counter',filename,timedepth_temp)
756	CALL Write_Ncdf_var('time_counter','time_counter', &
757	Child_file,timedepth_temp,'float')
758	CALL Copy_Ncdf_att('time_counter',TRIM(filename),Child_file)
759	DEALLOCATE(timedepth_temp)
760	varname = TRIM(Ncdf_varname(i))
761	Interpolation = .FALSE.
762	!
763	!copy time from input file to output file
764	CASE('time')
765	ALLOCATE(timedepth_temp(time))
766	CALL Read_Ncdf_var('time',filename,timedepth_temp)
767	CALL Write_Ncdf_var('time','time', &
768	Child_file,timedepth_temp,'float')
769	CALL Copy_Ncdf_att('time',TRIM(filename),Child_file)
770	DEALLOCATE(timedepth_temp)
771	varname = TRIM(Ncdf_varname(i))
772	Interpolation = .FALSE.
773	!
774	!copy deptht from input file to output file
775	CASE('deptht')
776	ALLOCATE(depth(deptht))
777	CALL Read_Ncdf_var('deptht',filename,depth)
778	CALL Write_Ncdf_var('deptht','deptht',Child_file,depth,'float')
779	CALL Copy_Ncdf_att('deptht',TRIM(filename),Child_file)
780	varname = TRIM(Ncdf_varname(i))
781	Interpolation = .FALSE.
782	!
783	!copy depthu from input file to output file
784	CASE('depthu')
785	ALLOCATE(depth(deptht))
786	CALL Read_Ncdf_var('depthu',filename,depth)
787	CALL Write_Ncdf_var('depthu','depthu',Child_file,depth,'float')
788	CALL Copy_Ncdf_att('depthu',TRIM(filename),Child_file)
789	varname = TRIM(Ncdf_varname(i))
790	Interpolation = .FALSE.
791	!
792	!copy depthv from input file to output file
793	CASE('depthv')
794	ALLOCATE(depth(deptht))
795	CALL Read_Ncdf_var('depthv',filename,depth)
796	CALL Write_Ncdf_var('depthv','depthv',Child_file,depth,'float')
797	CALL Copy_Ncdf_att('depthv',TRIM(filename),Child_file)
798	varname = TRIM(Ncdf_varname(i))
799	Interpolation = .FALSE.
800	!
801	!copy depthv from input file to output file
802	CASE('depthw')
803	ALLOCATE(depth(deptht))
804	CALL Read_Ncdf_var('depthw',filename,depth)
805	CALL Write_Ncdf_var('depthw','depthw',Child_file,depth,'float')
806	CALL Copy_Ncdf_att('depthw',TRIM(filename),Child_file)
807	varname = TRIM(Ncdf_varname(i))
808	Interpolation = .FALSE.
809	!
810	!copy time_steps from input file in case of restget use in NEMO in/out routines
811	CASE('time_steps')
812	! print *,'avant time step'
813
814	CALL Read_Ncdf_var('time_steps',filename,tabtemp1DInt)
815	! print *,'timedeph = ',tabtemp1DInt
816	CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt,'integer')
817	CALL Copy_Ncdf_att('time_steps',filename,Child_file)
818	DEALLOCATE(tabtemp1DInt)
819	Interpolation = .FALSE.
820	!
821	!store tmask in output file
822	CASE('tmask')
823	dimnames(1)='x'
824	dimnames(2)='y'
825	dimnames(3)='deptht'
826	IF (.NOT.ASSOCIATED(G1%tmask)) THEN
827	ALLOCATE(G1%tmask(nxfin,nyfin,deptht))
828	G1%tmask = 1
829	ENDIF
830	CALL Write_Ncdf_var('tmask',dimnames(1:3),Child_file,G1%tmask,'float')
831	varname = TRIM(Ncdf_varname(i))
832	Interpolation = .FALSE.
833	!
834	CASE default
835	varname = Ncdf_varname(i)
836	conservation = .FALSE.
837	CALL get_interptype( varname,interp_type,conservation )
838	WRITE(,) '**********************************************'
839	WRITE(,) 'varname = ',TRIM(varname), 'at ', posvar, ' point'
840	WRITE(,) 'interp_type = ',TRIM(interp_type)
841	WRITE(,) 'conservation = ',conservation
842	WRITE(,) '***********************************************'
843	Interpolation = .TRUE.
844	!
845	END SELECT
846
847	! //////////////// INTERPOLATION FOR 2D VARIABLES /////////////////////////////////////
848	!
849	IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 2 ) THEN
850	!
851	ALLOCATE(detected_pts(numlon,numlat,N))
852	ALLOCATE(masksrc(numlon,numlat))
853	!
854	IF(extrapolation) THEN
855	WRITE(,) 'interpolation/extrapolation ',TRIM(varname)
856	ELSE
857	WRITE(,) 'interpolation ',TRIM(varname)
858	ENDIF
859	!
860	ALLOCATE(tabvar3d(numlon,numlat,1))
861	ALLOCATE(tabinterp3d(nxfin,nyfin,1))
862	!
863	CALL Read_Ncdf_var(varname,filename,tabvar3d)
864	!
865	! search points where extrapolation is required
866	!
867	IF(Extrapolation) THEN
868	WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
869	CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
870	CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
871	ELSE
872	masksrc = .TRUE.
873	ENDIF
874
875	SELECT CASE(TRIM(interp_type))
876	CASE('bilinear')
877	CALL get_remap_matrix(latParent,latChild, &
878	lonParent,lonChild, &
879	masksrc,matrix,src_add,dst_add)
880
881	CASE('bicubic')
882	CALL get_remap_bicub(latParent,latChild, &
883	lonParent,lonChild, &
884	masksrc,matrix,src_add,dst_add)
885
886	END SELECT
887	!
888	!
889	SELECT CASE(TRIM(interp_type))
890	CASE('bilinear')
891	CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
892	matrix,src_add,dst_add)
893	CASE('bicubic')
894	CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
895	matrix,src_add,dst_add)
896	END SELECT
897	!
898	IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
899	G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)
900	!
901	IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)
902	!
903	dimnames(1)='x'
904	dimnames(2)='y'
905	!
906	CALL Write_Ncdf_var(TRIM(varname),dimnames(1:2),&
907	Child_file,tabinterp3d(:,:,1),'float')
908	!
909	DEALLOCATE(tabinterp3d)
910	DEALLOCATE(tabvar3d)
911	!
912	DEALLOCATE(detected_pts)
913	IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
914	DEALLOCATE( masksrc)
915
916	CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)
917	!
918
919	! //////////////// INTERPOLATION FOR 3D VARIABLES /////////////////////////////////////
920	!
921	ELSEIF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3 ) THEN
922	!
923	IF( DimUnlimited_Var(TRIM(varname),TRIM(filename)) ) THEN
924
925	ALLOCATE(detected_pts(numlon,numlat,N))
926	ALLOCATE(masksrc(numlon,numlat))
927	!
928	! ****************************LOOP ON TIME*****************************************
929	!loop on time
930	DO t = 1,time
931	!
932	IF(extrapolation) THEN
933	WRITE(,) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t
934	ELSE
935	WRITE(,) 'interpolation ',TRIM(varname),' for time t = ',t
936	ENDIF
937	!
938	ALLOCATE(tabvar3d(numlon,numlat,1))
939	ALLOCATE(tabinterp3d(nxfin,nyfin,1))
940	!
941	CALL Read_Ncdf_var(varname,filename,tabvar3d,t)
942	!
943	! search points where extrapolation is required
944	!
945	IF(Extrapolation) THEN
946	WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
947	IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
948	CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
949	ELSE
950	masksrc = .TRUE.
951	ENDIF
952
953	IF (t.EQ.1 ) THEN
954
955	SELECT CASE(TRIM(interp_type))
956	CASE('bilinear')
957	CALL get_remap_matrix(latParent,latChild, &
958	lonParent,lonChild, &
959	masksrc,matrix,src_add,dst_add)
960
961	CASE('bicubic')
962	CALL get_remap_bicub(latParent,latChild, &
963	lonParent,lonChild, &
964	masksrc,matrix,src_add,dst_add)
965
966	END SELECT
967	!
968	ENDIF
969	!
970	SELECT CASE(TRIM(interp_type))
971	CASE('bilinear')
972	CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
973	matrix,src_add,dst_add)
974	CASE('bicubic')
975	CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
976	matrix,src_add,dst_add)
977	END SELECT
978	!
979	IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
980	G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)
981	!
982	IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)
983	!
984	dimnames(1)='x'
985	dimnames(2)='y'
986	dimnames(3)='time_counter'
987	!
988	CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),&
989	Child_file,tabinterp3d,t,'float')
990	!
991	DEALLOCATE(tabinterp3d)
992	DEALLOCATE(tabvar3d)
993	!end loop on time
994	END DO
995	!
996	DEALLOCATE(detected_pts)
997	IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
998	DEALLOCATE( masksrc)
999
1000	ELSE
1001
1002	dimnames(1)='x'
1003	dimnames(2)='y'
1004	IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht'
1005	IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu'
1006	IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv'
1007	IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw'
1008	IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z'
1009	!
1010	! loop on vertical levels
1011	!
1012	DO nb = 1,deptht
1013	!
1014	ALLOCATE(masksrc(numlon,numlat))
1015	ALLOCATE(detected_pts(numlon,numlat,N))
1016	!
1017	! point detection et level n
1018	!
1019	land_level = .FALSE.
1020	IF( Extrapolation ) THEN
1021	IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE.
1022	ENDIF
1023
1024
1025	IF ( land_level ) THEN
1026	!
1027	WRITE(,) 'only land points on level ',nb
1028	ALLOCATE(tabinterp3d(nxfin,nyfin,1))
1029	tabinterp3d = 0.e0
1030	!
1031	CALL Write_Ncdf_var(TRIM(varname),dimnames, &
1032	Child_file,tabinterp3d,nb,'float')
1033	DEALLOCATE(tabinterp3d)
1034	!
1035	ELSE
1036	!
1037	ALLOCATE(tabvar01(numlon,numlat,1)) ! level k
1038	IF(Extrapolation) ALLOCATE(tabvar02(numlon,numlat,1)) ! level k-1
1039	IF(Extrapolation) ALLOCATE(tabvar03(numlon,numlat,1)) ! level k-2
1040	ALLOCATE(tabinterp3d(nxfin,nyfin,1))
1041	!
1042	IF(Extrapolation) THEN
1043	IF(nb==1) THEN
1044	CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
1045	WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
1046	ELSE IF (nb==2) THEN
1047	CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1)
1048	CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
1049	WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
1050	WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0.
1051	ELSE
1052	CALL Read_Ncdf_var(varname,filename,tabvar03,nb-2)
1053	CALL Read_Ncdf_var(varname,filename,tabvar02,nb-1)
1054	CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
1055	WHERE( tabvar01 .GE. 1.e+20 ) tabvar01 = 0.
1056	WHERE( tabvar02 .GE. 1.e+20 ) tabvar02 = 0.
1057	WHERE( tabvar03 .GE. 1.e+20 ) tabvar03 = 0.
1058	ENDIF
1059	!
1060	CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb)
1061
1062	ALLOCATE(tabvar1(numlon,numlat,1,1)) ! level k
1063	ALLOCATE(tabvar2(numlon,numlat,1,1)) ! level k-1
1064	ALLOCATE(tabvar3(numlon,numlat,1,1)) ! level k-2
1065	tabvar1(:,:,1,1) = tabvar01(:,:,1)
1066	tabvar2(:,:,1,1) = tabvar02(:,:,1)
1067	tabvar3(:,:,1,1) = tabvar03(:,:,1)
1068	CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,&
1069	tabvar3,G0,depth,masksrc,nb,'posvar')
1070	tabvar01(:,:,1) = tabvar1(:,:,1,1)
1071	tabvar02(:,:,1) = tabvar2(:,:,1,1)
1072	tabvar03(:,:,1) = tabvar3(:,:,1,1)
1073	DEALLOCATE(tabvar1,tabvar2,tabvar3)
1074	DEALLOCATE(tabvar02,tabvar03)
1075
1076	ELSE
1077	CALL Read_Ncdf_var(varname,filename,tabvar01,nb)
1078	IF(MAXVAL(tabvar01(:,:,1))==0.) land_level = .TRUE.
1079	masksrc = .TRUE.
1080	ENDIF
1081
1082	IF( Extrapolation ) THEN
1083	WRITE(,) 'interpolation/extrapolation ',TRIM(varname),' for vertical level = ',nb
1084	ELSE
1085	WRITE(,) 'interpolation ',TRIM(varname),' for vertical level = ',nb
1086	ENDIF
1087	!
1088	SELECT CASE(TRIM(interp_type))
1089	CASE('bilinear')
1090	CALL get_remap_matrix(latParent,latChild, &
1091	lonParent,lonChild, &
1092	masksrc,matrix,src_add,dst_add)
1093
1094	CASE('bicubic')
1095	CALL get_remap_bicub(latParent,latChild, &
1096	lonParent,lonChild, &
1097	masksrc,matrix,src_add,dst_add)
1098	!
1099	END SELECT
1100	!
1101	!
1102	SELECT CASE(TRIM(interp_type))
1103	!
1104	CASE('bilinear')
1105	CALL make_remap(tabvar01(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
1106	matrix,src_add,dst_add)
1107	CASE('bicubic')
1108	CALL make_bicubic_remap(tabvar01(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
1109	matrix,src_add,dst_add)
1110	END SELECT
1111	!
1112	IF( conservation ) CALL Correctforconservation(tabvar01(:,:,1),tabinterp3d(:,:,1), &
1113	G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)
1114
1115	!
1116	ALLOCATE(tabinterp4d(nxfin,nyfin,1,1))
1117	tabinterp4d(:,:,1,1) = tabinterp3d(:,:,1)
1118	IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb)
1119	tabinterp3d(:,:,1) = tabinterp4d(:,:,1,1)
1120	DEALLOCATE(tabinterp4d)
1121	!
1122	IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,nb)
1123	!
1124	CALL Write_Ncdf_var(TRIM(varname),dimnames, &
1125	Child_file,tabinterp3d,nb,'float')
1126	!
1127	DEALLOCATE(tabinterp3d)
1128	DEALLOCATE(tabvar01)
1129	!
1130	!
1131	ENDIF
1132
1133	!
1134	IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
1135	DEALLOCATE( masksrc )
1136	DEALLOCATE(detected_pts)
1137	!
1138	! end loop on vertical levels
1139	!
1140	END DO
1141
1142
1143	ENDIF
1144
1145	CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)
1146	!
1147	ELSE IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 4 ) THEN
1148	!
1149	!
1150	! //////////////// INTERPOLATION FOR 4D VARIABLES /////////////////////////////////////
1151	!
1152	dimnames(1)='x'
1153	dimnames(2)='y'
1154	IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht'
1155	IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu'
1156	IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv'
1157	IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw'
1158	IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z'
1159	dimnames(4)='time_counter'
1160	!
1161	! loop on vertical levels
1162	!
1163	DO nb = 1,deptht
1164	!
1165	ALLOCATE(masksrc(numlon,numlat))
1166	ALLOCATE(detected_pts(numlon,numlat,N))
1167	!
1168	! point detection et level n
1169	!
1170	land_level = .FALSE.
1171	IF( Extrapolation ) THEN
1172	IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE.
1173	ENDIF
1174
1175
1176	IF ( land_level ) THEN
1177	!
1178	WRITE(,) 'only land points on level ',nb
1179	ALLOCATE(tabinterp4d(nxfin,nyfin,1,1))
1180	tabinterp4d = 0.e0
1181	!
1182	DO ii = 1,time
1183	CALL Write_Ncdf_var(TRIM(varname),dimnames, &
1184	Child_file,tabinterp4d,ii,nb,'float')
1185	END DO
1186	DEALLOCATE(tabinterp4d)
1187	!
1188	ELSE
1189	!
1190	! loop on time
1191	!
1192	DO t = 1,time
1193	!
1194	ALLOCATE(tabvar1(numlon,numlat,1,1)) ! level k
1195	IF(Extrapolation) ALLOCATE(tabvar2(numlon,numlat,1,1)) ! level k-1
1196	IF(Extrapolation) ALLOCATE(tabvar3(numlon,numlat,1,1)) ! level k-2
1197	ALLOCATE(tabinterp4d(nxfin,nyfin,1,1))
1198	!
1199	IF(Extrapolation) THEN
1200	IF(nb==1) THEN
1201	CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
1202	WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
1203	ELSE IF (nb==2) THEN
1204	CALL Read_Ncdf_var(varname,filename,tabvar2,t,nb-1)
1205	CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
1206	WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
1207	WHERE( tabvar2 .GE. 1.e+20 ) tabvar2 = 0.
1208	ELSE
1209	CALL Read_Ncdf_var(varname,filename,tabvar3,t,nb-2)
1210	CALL Read_Ncdf_var(varname,filename,tabvar2,t,nb-1)
1211	CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
1212	WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
1213	WHERE( tabvar2 .GE. 1.e+20 ) tabvar2 = 0.
1214	WHERE( tabvar3 .GE. 1.e+20 ) tabvar3 = 0.
1215	ENDIF
1216	!
1217	IF (t.EQ.1 ) CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb)
1218
1219	CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,&
1220	tabvar3,G0,depth,masksrc,nb,'posvar')
1221	DEALLOCATE(tabvar2,tabvar3)
1222
1223	ELSE
1224	CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
1225	IF(MAXVAL(tabvar1(:,:,1,1))==0.) land_level = .TRUE.
1226	masksrc = .TRUE.
1227	ENDIF
1228
1229	IF( Extrapolation ) THEN
1230	WRITE(,) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb
1231	ELSE
1232	WRITE(,) 'interpolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb
1233	ENDIF
1234
1235	!
1236
1237	IF (t.EQ.1 ) THEN
1238
1239	SELECT CASE(TRIM(interp_type))
1240	CASE('bilinear')
1241	CALL get_remap_matrix(latParent,latChild, &
1242	lonParent,lonChild, &
1243	masksrc,matrix,src_add,dst_add)
1244
1245	CASE('bicubic')
1246	CALL get_remap_bicub(latParent,latChild, &
1247	lonParent,lonChild, &
1248	masksrc,matrix,src_add,dst_add)
1249	!
1250	END SELECT
1251	!
1252	ENDIF
1253	!
1254	SELECT CASE(TRIM(interp_type))
1255	!
1256	CASE('bilinear')
1257	CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
1258	matrix,src_add,dst_add)
1259	CASE('bicubic')
1260	CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),nxfin,nyfin, &
1261	matrix,src_add,dst_add)
1262	END SELECT
1263	!
1264	IF( conservation ) CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
1265	G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin)
1266
1267	!
1268	IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb)
1269	!
1270	IF(Extrapolation) tabinterp4d(:,:,1,1) = tabinterp4d(:,:,1,1) * mask(:,:,nb)
1271	!
1272	CALL Write_Ncdf_var(TRIM(varname),dimnames, &
1273	Child_file,tabinterp4d,t,nb,'float')
1274	!
1275	DEALLOCATE(tabinterp4d)
1276	DEALLOCATE(tabvar1)
1277	!
1278	! end loop on time
1279	!
1280
1281	END DO
1282
1283	ENDIF
1284
1285	!
1286	IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
1287	DEALLOCATE( masksrc )
1288	DEALLOCATE(detected_pts)
1289	!
1290	! end loop on vertical levels
1291	!
1292	END DO
1293	!
1294	CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)
1295	!
1296	! fin du if interpolation ...
1297	!
1298	ENDIF
1299	!
1300	END DO
1301
1302	PRINT *,'FIN DE INTERPEXTRAPVAR'
1303	!
1304	IF(Extrapolation) DEALLOCATE(G1%tmask,G0%tmask)
1305	DEALLOCATE(G0%e1t,G0%e2t,G1%e1t,G1%e2t)
1306	IF(ASSOCIATED(depth)) DEALLOCATE(depth)
1307	!
1308	END SUBROUTINE Interp_Extrap_var
1309	!
1310	!
1311	END MODULE agrif_interpolation

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_r4.0-HEAD/NESTING/src/agrif_interpolation.f90 @ 14974

Download in other formats: