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

Last change on this file since 12253 was 12243, checked in by clem, 12 months ago

avoid out of bounds when the size of the regional domain is exactly the same as the size of the reference domain (with ln_agrif_domain=F)

  • Property svn:keywords set to Id
File size: 36.5 KB
Line 
1!************************************************************************
2! Fortran 95 OPA Nesting tools                  *
3!                          *
4!     Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr)   *
5!                          *
6!************************************************************************
7!
8MODULE agrif_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  !
29CONTAINS
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( zx*jif-1, zx*irafx ) / REAL(zx*irafx)
249          By = MOD( zy*jjf-1, zy*irafy ) / 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 = k*s(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 => 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
519
520    !       
521    TYPE(Coordinates) :: G0,G1     
522    !
523    !*****************
524    !If coarse grid is masked possibility to activate an extrapolation process
525    !
526    Extrapolation = .FALSE.
527    PRINT*,'DEBUT INTERP_EXTRAP_VAR'
528    !
529    !*****************
530    !
531    ! check grid position
532    !
533    IF( PRESENT(cd_type) ) THEN
534       posvar = cd_type
535    ELSE
536       posvar = 'T'
537    ENDIF     
538
539    !
540    ! read dimensions in netcdf file
541    !
542    CALL Read_Ncdf_dim('x',filename,numlon)
543    CALL Read_Ncdf_dim('y',filename,numlat)
544    CALL Read_Ncdf_dim('time_counter',filename,time)
545    IF ( Dims_Existence( 'deptht' , filename ) ) THEN
546       CALL Read_Ncdf_dim('deptht',filename,deptht)
547    ELSE IF ( Dims_Existence( 'depthu' , filename ) ) THEN
548       CALL Read_Ncdf_dim('depthu',filename,deptht)
549    ELSE IF ( Dims_Existence( 'depthv' , filename ) ) THEN
550       CALL Read_Ncdf_dim('depthv',filename,deptht)
551    ELSE IF ( Dims_Existence( 'depthw' , filename ) ) THEN
552       CALL Read_Ncdf_dim('depthw',filename,deptht)
553    ELSE IF ( Dims_Existence( 'z' , filename ) ) THEN
554       CALL Read_Ncdf_dim('z',filename,deptht)
555    ELSE
556       deptht = N
557    ENDIF
558
559    !
560    ! retrieve netcdf variable name
561    !
562    CALL Read_Ncdf_VarName(filename,Ncdf_varname)
563    !
564    ! define name of child grid file
565    !
566    CALL set_child_name(filename,Child_file) 
567    CALL set_child_name(parent_coordinate_file,Childcoordinates)
568    CALL set_child_name(parent_bathy_level,childbathy)
569    WRITE(*,*) 'Child grid file name = ',TRIM(Child_file)
570    !
571    ! create this file
572    !
573    status = nf90_create(Child_file,NF90_WRITE,ncid)
574    status = nf90_close(ncid)
575    !
576    ! read coordinates of both domains
577    !           
578    status = Read_Coordinates(parent_coordinate_file,G0)
579    status = Read_Coordinates(Childcoordinates,G1,Pacifique)
580    !
581    ! check consistency of informations read in namelist
582    !     
583    IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
584         imax <= imin .OR. jmax <= jmin ) THEN                   
585       WRITE(*,*) 'ERROR ***** bad child grid definition ...'
586       WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
587       WRITE(*,*) ' '
588       STOP
589    ENDIF
590    !
591    IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
592       WRITE(*,*) 'ERROR ***** bad child coordinates or bathy file ...'
593       WRITE(*,*) ' '
594       WRITE(*,*) 'please check that child coordinates file and child bathymetry file'
595       WRITE(*,*) 'has been created with the current namelist '
596       WRITE(*,*) ' '
597       STOP
598    ENDIF
599    !     
600
601    !
602    ! Initialization of T-mask thanks to bathymetry
603    !
604    IF( Extrapolation ) THEN
605
606       WRITE(*,*) 'mask initialisation on coarse and fine grids'
607       !
608       ALLOCATE(mask(numlon,numlat,N))
609       CALL Init_mask(childbathy,G1,1,1)
610       CALL Init_mask(parent_bathy_level,G0,1,1)
611       !
612    ENDIF
613
614    !     
615    ! select coordinates to use according to variable position
616    !
617    ALLOCATE(lonParent(numlon,numlat),latParent(numlon,numlat))
618    ALLOCATE(lonChild(nxfin,nyfin),latChild(nxfin,nyfin))
619
620    SELECT CASE(posvar)
621    CASE('T')
622       lonParent = G0%glamt
623       latParent = G0%gphit
624       lonChild  = G1%glamt
625       latChild  = G1%gphit
626       mask      = G1%tmask
627    CASE('U')
628       lonParent = G0%glamu
629       latParent = G0%gphiu
630       lonChild  = G1%glamu
631       latChild  = G1%gphiu
632       mask      = G1%umask
633    CASE('V')
634       lonParent = G0%glamv
635       latParent = G0%gphiv
636       lonChild  = G1%glamv
637       latChild  = G1%gphiv
638       mask      = G1%vmask
639    END SELECT
640
641    DEALLOCATE(G0%glamu,G0%glamv,G0%gphiu,G0%gphiv)
642    DEALLOCATE(G1%glamu,G1%glamv,G1%gphiu,G1%gphiv) 
643    DEALLOCATE(G1%glamt,G1%gphit,G0%glamt,G0%gphit)
644
645    !
646    !longitude modification if child domain covers Pacific ocean area
647    !
648    IF( lonChild(1,1) > lonChild(nxfin,nyfin) ) THEN
649       Pacifique = .TRUE.
650       WHERE( lonChild < 0 )
651          lonChild = lonChild + 360.
652       END WHERE
653       WHERE( lonParent < 0 )
654          lonParent = lonParent + 360.
655       END WHERE
656    ENDIF
657
658    !   
659    !
660    ! dimensions initialization     
661    !
662    CALL Write_Ncdf_dim('x',Child_file,nxfin)
663    CALL Write_Ncdf_dim('y',Child_file,nyfin)
664    IF ( Dims_Existence( 'deptht' , filename ) ) CALL Write_Ncdf_dim('deptht',Child_file,deptht)
665    IF ( Dims_Existence( 'depthu' , filename ) ) CALL Write_Ncdf_dim('depthu',Child_file,deptht)
666    IF ( Dims_Existence( 'depthv' , filename ) ) CALL Write_Ncdf_dim('depthv',Child_file,deptht)
667    IF ( Dims_Existence( 'depthw' , filename ) ) CALL Write_Ncdf_dim('depthw',Child_file,deptht)
668    IF ( Dims_Existence( 'z' , filename ) ) CALL Write_Ncdf_dim('z',Child_file,deptht)
669    CALL Write_Ncdf_dim('time_counter',Child_file,0)
670
671    IF( deptht .NE. 1 .AND. deptht .NE. N ) THEN
672       WRITE(*,*) '***'
673       WRITE(*,*) 'Number of vertical levels doesn t match between namelist'
674       WRITE(*,*) 'and forcing file ',TRIM(filename)
675       WRITE(*,*) 'Please check the values in namelist file'
676       WRITE(*,*) 'N = ',N
677       WRITE(*,*) 'deptht = ',deptht
678       STOP     
679    ENDIF
680    !
681    !
682    DO i = 1,SIZE(Ncdf_varname)
683       !     
684       !
685       ! ******************************LOOP ON VARIABLE NAMES*******************************************
686       !
687       !     
688       SELECT CASE (TRIM(Ncdf_varname(i)))
689          !
690          !copy nav_lon from child coordinates to output file     
691       CASE('nav_lon')   
692          ALLOCATE(latlon_temp(nxfin,nyfin))
693          CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),latlon_temp) 
694          CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,latlon_temp,'float')
695          CALL Copy_Ncdf_att('nav_lon',TRIM(filename),Child_file, &
696               MINVAL(latlon_temp),MAXVAL(latlon_temp))
697          DEALLOCATE(latlon_temp)
698          varname = TRIM(Ncdf_varname(i))
699          Interpolation = .FALSE.
700          !     
701          !copy nav_lat from child coordinates to output file
702       CASE('nav_lat')             
703          ALLOCATE(latlon_temp(nxfin,nyfin))
704          CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),latlon_temp) 
705          CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,latlon_temp,'float')
706          CALL Copy_Ncdf_att('nav_lat',TRIM(filename),Child_file, &
707               MINVAL(latlon_temp),MAXVAL(latlon_temp)) 
708          DEALLOCATE(latlon_temp)
709          varname = TRIM(Ncdf_varname(i))
710          Interpolation = .FALSE.
711          !
712          !copy nav_lev from restart_file to output file
713          !
714       CASE('nav_lev')
715
716          WRITE(*,*) 'copy nav_lev'
717          CALL Read_Ncdf_var('nav_lev',filename,nav_lev) 
718          IF(.NOT. dimg ) THEN
719             CALL Write_Ncdf_var('nav_lev','z',Child_file,nav_lev,'float')
720             CALL Copy_Ncdf_att('nav_lev',filename,Child_file)     
721          ENDIF
722          DEALLOCATE(nav_lev)
723          Interpolation = .FALSE.
724          !       
725          !copy time_counter from input file to output file
726       CASE('time_counter')
727          ALLOCATE(timedepth_temp(time))
728          CALL Read_Ncdf_var('time_counter',filename,timedepth_temp) 
729          CALL Write_Ncdf_var('time_counter','time_counter',  &
730               Child_file,timedepth_temp,'float')
731          CALL Copy_Ncdf_att('time_counter',TRIM(filename),Child_file)
732          DEALLOCATE(timedepth_temp)
733          varname = TRIM(Ncdf_varname(i))
734          Interpolation = .FALSE.
735          !
736          !copy deptht from input file to output file
737       CASE('deptht')
738          ALLOCATE(depth(deptht))
739          CALL Read_Ncdf_var('deptht',filename,depth) 
740          CALL Write_Ncdf_var('deptht','deptht',Child_file,depth,'float')
741          CALL Copy_Ncdf_att('deptht',TRIM(filename),Child_file)
742          varname = TRIM(Ncdf_varname(i))
743          Interpolation = .FALSE.     
744          !
745          !copy depthu from input file to output file
746       CASE('depthu')
747          ALLOCATE(depth(deptht))
748          CALL Read_Ncdf_var('depthu',filename,depth) 
749          CALL Write_Ncdf_var('depthu','depthu',Child_file,depth,'float')
750          CALL Copy_Ncdf_att('depthu',TRIM(filename),Child_file)
751          varname = TRIM(Ncdf_varname(i))
752          Interpolation = .FALSE.     
753          !
754          !copy depthv from input file to output file
755       CASE('depthv')
756          ALLOCATE(depth(deptht))
757          CALL Read_Ncdf_var('depthv',filename,depth) 
758          CALL Write_Ncdf_var('depthv','depthv',Child_file,depth,'float')
759          CALL Copy_Ncdf_att('depthv',TRIM(filename),Child_file)
760          varname = TRIM(Ncdf_varname(i))
761          Interpolation = .FALSE.     
762          !
763          !copy depthv from input file to output file
764       CASE('depthw')
765          ALLOCATE(depth(deptht))
766          CALL Read_Ncdf_var('depthw',filename,depth) 
767          CALL Write_Ncdf_var('depthw','depthw',Child_file,depth,'float')
768          CALL Copy_Ncdf_att('depthw',TRIM(filename),Child_file)
769          varname = TRIM(Ncdf_varname(i))
770          Interpolation = .FALSE.     
771          !
772          !copy time_steps from input file in case of restget use in NEMO in/out routines
773       CASE('time_steps')
774          !         print *,'avant time step'
775
776          CALL Read_Ncdf_var('time_steps',filename,tabtemp1DInt) 
777          !       print *,'timedeph = ',tabtemp1DInt
778          CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt,'integer')
779          CALL Copy_Ncdf_att('time_steps',filename,Child_file) 
780          DEALLOCATE(tabtemp1DInt)
781          Interpolation = .FALSE. 
782         
783          !store tmask in output file
784       CASE('tmask')
785          dimnames(1)='x'
786          dimnames(2)='y'
787          dimnames(3)='deptht' 
788          IF (.NOT.ASSOCIATED(G1%tmask)) THEN
789             ALLOCATE(G1%tmask(nxfin,nyfin,deptht))
790             G1%tmask = 1
791          ENDIF
792          CALL Write_Ncdf_var('tmask',dimnames(1:3),Child_file,G1%tmask,'float')
793          varname = TRIM(Ncdf_varname(i))
794          Interpolation = .FALSE.     
795          !     
796       CASE default
797          varname = Ncdf_varname(i)
798          conservation = .FALSE.
799          CALL get_interptype( varname,interp_type,conservation )
800          WRITE(*,*) '**********************************************'
801          WRITE(*,*) 'varname      = ',TRIM(varname), 'at ', posvar, ' point'
802          WRITE(*,*) 'interp_type  = ',TRIM(interp_type)
803          WRITE(*,*) 'conservation = ',conservation
804          WRITE(*,*) '***********************************************'
805          Interpolation = .TRUE.
806          !         
807       END SELECT
808
809       ! //////////////// INTERPOLATION FOR 3D VARIABLES /////////////////////////////////////
810       !
811       IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3 ) THEN
812          !     
813          ALLOCATE(detected_pts(numlon,numlat,N))
814          ALLOCATE(masksrc(numlon,numlat))                                           
815          !
816          ! ******************************LOOP ON TIME*******************************************
817          !loop on time
818          DO t = 1,time
819             !                   
820             IF(extrapolation) THEN
821                WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t
822             ELSE
823                WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t
824             ENDIF
825             !                           
826             ALLOCATE(tabvar3d(numlon,numlat,1))     
827             ALLOCATE(tabinterp3d(nxfin,nyfin,1))
828             !                   
829             CALL Read_Ncdf_var(varname,filename,tabvar3d,t)                                             
830             !
831             ! search points where extrapolation is required
832             !
833             IF(Extrapolation) THEN
834                 WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
835                IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
836                CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
837             ELSE
838                masksrc = .TRUE.
839             ENDIF
840
841             IF (t.EQ.1 ) THEN
842
843                SELECT CASE(TRIM(interp_type))
844                CASE('bilinear')
845                   CALL get_remap_matrix(latParent,latChild,   &
846                        lonParent,lonChild,                    &
847                        masksrc,matrix,src_add,dst_add)
848
849                CASE('bicubic')
850                   CALL get_remap_bicub(latParent,latChild,   &
851                        lonParent,lonChild,   &
852                        masksrc,matrix,src_add,dst_add)
853
854                END SELECT
855                !                                                       
856             ENDIF
857             !     
858             SELECT CASE(TRIM(interp_type))
859             CASE('bilinear')                                                       
860                CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
861                     matrix,src_add,dst_add)     
862             CASE('bicubic')                                   
863                CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
864                     matrix,src_add,dst_add)                       
865             END SELECT
866             !                     
867             IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
868                  G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 
869             !     
870             IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)     
871             !                     
872             dimnames(1)='x'
873             dimnames(2)='y'
874             dimnames(3)='time_counter'
875             !                                     
876             CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),&
877                  Child_file,tabinterp3d,t,'float')
878             !     
879             DEALLOCATE(tabinterp3d)
880             DEALLOCATE(tabvar3d)                     
881             !end loop on time
882          END DO
883          !                   
884          DEALLOCATE(detected_pts)
885          IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
886          DEALLOCATE( masksrc)
887
888          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)                       
889          !
890       ELSE IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 4 ) THEN
891          !
892          !
893          ! //////////////// INTERPOLATION FOR 4D VARIABLES /////////////////////////////////////
894          !
895          dimnames(1)='x'
896          dimnames(2)='y'
897          IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht' 
898          IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu'
899          IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv'
900          IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw'
901          IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z'
902          dimnames(4)='time_counter' 
903          !
904          ! loop on vertical levels
905          !
906          DO nb = 1,deptht
907             !       
908             ALLOCATE(masksrc(numlon,numlat))
909             ALLOCATE(detected_pts(numlon,numlat,N))         
910             !
911             ! point detection et level n                                     
912             !
913             land_level = .FALSE.
914             IF( Extrapolation ) THEN
915                IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE.
916             ENDIF
917
918
919             IF ( land_level ) THEN
920                !
921                WRITE(*,*) 'only land points on level ',nb
922                ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
923                tabinterp4d = 0.e0
924                !                             
925                DO ii = 1,time
926                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
927                        Child_file,tabinterp4d,ii,nb,'float')
928                END DO
929                DEALLOCATE(tabinterp4d)
930                !
931             ELSE
932                !                       
933                ! loop on time                       
934                !
935                DO t = 1,time
936                   !                   
937                   ALLOCATE(tabvar1(numlon,numlat,1,1))    ! level k
938                   IF(Extrapolation) ALLOCATE(tabvar2(numlon,numlat,1,1))    ! level k-1
939                   IF(Extrapolation) ALLOCATE(tabvar3(numlon,numlat,1,1))    ! level k-2                       
940                   ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
941                   !
942                   IF(Extrapolation) THEN                                 
943                      IF(nb==1) THEN
944                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)   
945                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
946                      ELSE IF (nb==2) THEN
947                         CALL Read_Ncdf_var(varname,filename,tabvar2,t,nb-1)
948                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)           
949                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
950                         WHERE( tabvar2 .GE. 1.e+20 ) tabvar2 = 0.
951                      ELSE
952                         CALL Read_Ncdf_var(varname,filename,tabvar3,t,nb-2)
953                         CALL Read_Ncdf_var(varname,filename,tabvar2,t,nb-1)
954                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
955                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
956                         WHERE( tabvar2 .GE. 1.e+20 ) tabvar2 = 0.
957                         WHERE( tabvar3 .GE. 1.e+20 ) tabvar3 = 0.
958                      ENDIF
959                      !                             
960                      IF (t.EQ.1 ) CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb)
961
962                      CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,&
963                           tabvar3,G0,depth,masksrc,nb,'posvar')                       
964                      DEALLOCATE(tabvar2,tabvar3)
965
966                   ELSE
967                      CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
968                      IF(MAXVAL(tabvar1(:,:,1,1))==0.) land_level = .TRUE.
969                      masksrc = .TRUE. 
970                   ENDIF
971
972                   IF( Extrapolation ) THEN
973                      WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb   
974                   ELSE
975                      WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb
976                   ENDIF
977
978                   !                         
979
980                   IF (t.EQ.1 ) THEN
981
982                      SELECT CASE(TRIM(interp_type))
983                      CASE('bilinear')
984                         CALL get_remap_matrix(latParent,latChild,   &
985                              lonParent,lonChild,   &
986                              masksrc,matrix,src_add,dst_add)
987
988                      CASE('bicubic')
989                         CALL get_remap_bicub(latParent,latChild,   &
990                              lonParent,lonChild,   &
991                              masksrc,matrix,src_add,dst_add)
992                         !                             
993                      END SELECT
994                      !                                                       
995                   ENDIF
996                   !     
997                   SELECT CASE(TRIM(interp_type))
998                      !                             
999                   CASE('bilinear')                                                       
1000                      CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
1001                           matrix,src_add,dst_add)     
1002                   CASE('bicubic')                                   
1003                      CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),nxfin,nyfin, &
1004                           matrix,src_add,dst_add)                       
1005                   END SELECT
1006                   !                     
1007                   IF( conservation ) CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
1008                        G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 
1009
1010                   !
1011                   IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb)
1012                   !     
1013                   IF(Extrapolation) tabinterp4d(:,:,1,1) = tabinterp4d(:,:,1,1) * mask(:,:,nb)     
1014                   !                     
1015                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
1016                        Child_file,tabinterp4d,t,nb,'float')
1017                   !     
1018                   DEALLOCATE(tabinterp4d)
1019                   DEALLOCATE(tabvar1)                     
1020                   !
1021                   ! end loop on time
1022                   !
1023
1024                END DO
1025
1026             ENDIF
1027
1028             !
1029             IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
1030             DEALLOCATE( masksrc )
1031             DEALLOCATE(detected_pts)       
1032             !
1033             ! end loop on vertical levels
1034             !     
1035          END DO
1036          !   
1037          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)
1038          !
1039          ! fin du if interpolation ...
1040          !
1041       ENDIF
1042       !
1043    END DO
1044
1045    PRINT *,'FIN DE INTERPEXTRAPVAR'
1046    !   
1047    IF(Extrapolation) DEALLOCATE(G1%tmask,G0%tmask)
1048    DEALLOCATE(G0%e1t,G0%e2t,G1%e1t,G1%e2t)
1049    IF(ASSOCIATED(depth)) DEALLOCATE(depth)   
1050    !
1051  END SUBROUTINE Interp_Extrap_var
1052  !
1053  !
1054END MODULE agrif_interpolation
Note: See TracBrowser for help on using the repository browser.