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_interpolation.f90 in utils/tools/NESTING/src – NEMO

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

Last change on this file was 12243, checked in by clem, 4 years 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.