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 @ 10025

Last change on this file since 10025 was 10025, checked in by clem, 6 years ago

nesting tools are partly rewritten (mostly for create_coordinates and bathy) to get better functionality. Now you can use the nesting to either define an agrif zoom or a regional domain (for bdy purposes). Also, the nesting tools output a domain_cfg.nc that can be directly used in NEMO4 without the need of DOMAINcfg tool. The option of median average for bathymetry interpolation still does not work properly but it's not new

  • Property svn:keywords set to Id
File size: 36.4 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,zx,zy
201    INTEGER :: ji,jj,jif,jjf,jic,jjc,jdecx,jdecy
202    REAL*8 :: Ax, Bx, Ay, By
203
204    nxf = SIZE(tabout,1)
205    nyf = SIZE(tabout,2)
206
207    SELECT CASE(typevar)
208    CASE('T')
209       IF(MOD(irafx,2)==1) THEN ! odd
210          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
211       ELSE                     ! even
212          zx = 2 ; zy = 2 ; jdecx = FLOOR(irafx/2.) ; jdecy = FLOOR(irafy/2.)
213       ENDIF
214    CASE('U')
215       IF(MOD(irafx,2)==1) THEN ! odd
216          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
217       ELSE                     ! even
218          zx = 1 ; zy = 2 ; jdecx = irafx - 1 ; jdecy = FLOOR(irafy/2.)
219       ENDIF
220    CASE('V')
221       IF(MOD(irafx,2)==1) THEN ! odd
222          zx = 1 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
223       ELSE                     ! even
224          zx = 2 ; zy = 1 ; jdecx = FLOOR(irafx/2.) ; jdecy = irafy - 1
225       ENDIF
226    CASE('F')
227       IF(MOD(irafx,2)==1) THEN ! odd
228          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
229       ELSE                     ! even
230          zx = 1 ; zy = 1 ; jdecx = irafx - 1 ; jdecy = irafy - 1
231       ENDIF
232    END SELECT
233
234   
235    DO jj = 1, nyf
236
237       jjf = jj - jdecy
238       jjc = j_min + FLOOR((jjf-1.) / irafy) 
239       
240       DO ji = 1, nxf
241         
242          jif = ji - jdecx 
243          jic = i_min + FLOOR((jif-1.) / irafx) 
244         
245          Bx = MOD( zx*jif-1, zx*irafx ) / REAL(zx*irafx)
246          By = MOD( zy*jjf-1, zy*irafy ) / REAL(zy*irafy)
247          Ax = 1. - Bx
248          Ay = 1. - By
249
250          tabout(ji,jj) = ( Bx * tabin(jic+1,jjc  ) + Ax * tabin(jic,jjc  ) ) * Ay + &
251             &            ( Bx * tabin(jic+1,jjc+1) + Ax * tabin(jic,jjc+1) ) * By
252       END DO
253    END DO
254   
255    !
256  END SUBROUTINE agrif_base_interp3
257
258  !       
259  SUBROUTINE polint(xin,valin,n,x,val)
260    IMPLICIT NONE
261    INTEGER n
262    REAL*8 xin(n), valin(n)
263    REAL*8 x,val
264
265    INTEGER ns,i,m
266    REAL *8 dif,dift
267    REAL*8 c(n),d(n),ho,hp,w,den,dy
268
269    ns = 1
270    dif = ABS(x-xin(1))
271
272    DO i=1,n
273       dift = ABS(x-xin(i))
274       IF (dift < dif) THEN
275     ns = i
276     dif = dift
277       ENDIF
278       c(i) = valin(i)
279       d(i) = valin(i)
280    ENDDO
281
282    val = valin(ns)
283    ns = ns - 1
284
285    DO m=1,n-1
286       DO i=1,n-m
287     ho = xin(i)-x
288     hp = xin(i+m)-x
289     w=c(i+1)-d(i)
290     den = w/(ho-hp)
291          d(i) = hp * den
292     c(i) = ho * den
293       ENDDO
294       IF (2*ns < (n-m)) THEN
295          dy = c(ns+1)
296       ELSE
297          dy = d(ns)
298          ns = ns - 1
299       ENDIF
300       val = val + dy
301    ENDDO
302  END SUBROUTINE polint
303  !
304  !     
305  SUBROUTINE polcoe(xin,valin,n,cof)
306    IMPLICIT NONE
307    INTEGER n
308    REAL*8 xin(n),valin(n),cof(n)
309
310
311    INTEGER i,j,k
312    REAL*8 b,ff,phi,s(n)
313
314    s = 0.
315    cof = 0.
316
317    s(n)=-xin(1)
318    DO i=2,n
319       DO j=n+1-i,n-1
320          s(j) =s(j) -xin(i)*s(j+1)
321       ENDDO
322       s(n)=s(n)-xin(i)
323    ENDDO
324
325    DO j=1,n
326       phi=n
327       DO k=n-1,1,-1
328          phi = k*s(k+1)+xin(j)*phi
329       ENDDO
330       ff=valin(j)/phi
331       b=1.
332       DO k=n,1,-1
333     cof(k)=cof(k)+b*ff
334     b=s(k)+xin(j)*b
335       ENDDO
336    ENDDO
337
338    RETURN
339  END SUBROUTINE polcoe
340  !
341
342  !****************************************************************
343  !   subroutine Correctforconservation            *
344  !                        *
345  ! Conservation on domain boundaries ( over 2 coarse grid cells) *
346  !                        *
347  !****************************************************************         
348  !
349  !
350  SUBROUTINE Correctforconservation(tabcoarse,tabfine,e1parent,e2parent,e1,e2,nxfin,nyfin,posvar,i_min,j_min)
351    IMPLICIT NONE
352    INTEGER :: nxfin,nyfin
353    REAL*8,DIMENSION(:,:) :: tabcoarse,tabfine,e1parent,e2parent,e1,e2
354    CHARACTER*1  :: posvar
355    INTEGER :: i_min,j_min,diff
356    INTEGER ji,jj,ipt,jpt,i,j
357    INTEGER ind1,ind2,ind3,ind4,ind5,ind6
358    REAL*8 cff1,cff2,cff3     
359    !
360    diff = 0     
361    IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
362    !       
363    ind1 = 2 + CEILING(irafx/2.0) + diff
364    ind2 = nxfin-(2-1)-CEILING(irafx/2.0)
365    ind3 = 2 + CEILING(irafy/2.0) + diff
366    ind4 = nyfin-(2-1)-CEILING(irafy/2.0)
367    ind5 = nxfin - 1 - irafx - CEILING(irafx/2.0) 
368    ind6 = nyfin - 1 - irafy - CEILING(irafy/2.0) 
369    !       
370    IF (posvar.EQ.'T') THEN
371       !       
372       PRINT *,'correction'
373       !
374       DO ji=ind1,ind1+irafx,irafx
375          !       
376          ipt=i_min+(3-1)+(ji-ind1)/irafx
377          !           
378          DO jj=ind3,ind4,irafy
379             !           
380             jpt=j_min+(3-1)+(jj-ind3)/irafy
381
382             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)* &
383                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
384
385             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)*   &
386                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
387                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
388
389             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                       
390
391             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
392                  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
393
394          END DO
395
396       ENDDO
397       !****     
398       DO ji=ind5,ind5+irafx,irafx
399          !       
400          ipt=i_min+(3-1)+(ji-ind1)/irafx
401
402          DO jj=ind3,ind4,irafy
403             jpt=j_min+(3-1)+(jj-ind3)/irafy
404
405             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)* &
406                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
407
408             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)*   &
409                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
410                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
411
412             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                       
413
414             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
415                  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
416
417          END DO
418       ENDDO
419
420       DO jj=ind3,ind3+irafy,irafy
421          !       
422          jpt=j_min+(3-1)+(jj-ind3)/irafy
423          !
424          DO ji=ind1,ind2,irafx
425             !         
426             ipt=i_min+(3-1)+(ji-ind1)/irafx
427
428             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)* &
429                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
430
431             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)*   &
432                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
433                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
434
435             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                       
436
437             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
438                  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
439
440          END DO
441       ENDDO
442       !****     
443       DO jj=ind6,ind6+irafy,irafy
444          !
445          jpt=j_min+(3-1)+(jj-ind3)/irafy
446          !
447          DO ji=ind1,ind2,irafx         
448             ipt=i_min+(3-1)+(ji-ind1)/irafx 
449
450             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)* &
451                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
452
453             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)*   &
454                  e2(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)*   &
455                  tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff))
456
457             cff3=e1parent(ipt,jpt)*e2parent(ipt,jpt)*tabcoarse(ipt,jpt)                       
458
459             tabfine(ji-FLOOR(irafx/2.0):ji+FLOOR(irafx/2.0)-diff,jj-FLOOR(irafy/2.0):jj+FLOOR(irafy/2.0)-diff)=  &
460                  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
461
462          END DO
463       ENDDO
464       !                                     
465       !       
466    ENDIF
467    RETURN
468  END SUBROUTINE Correctforconservation
469
470  !
471  !****************************************************************
472  !   subroutine Interp_Extrap_var           *
473  !                        *
474  ! Interpolation and Extrapolation for temperature and salinity  *
475  !                        *
476  ! -input :                     *
477  !       filename : file containing variable to interpolate   *
478  !                        *
479  !****************************************************************         
480  !
481  SUBROUTINE Interp_Extrap_var(filename, cd_type)     
482    !
483    IMPLICIT NONE
484    !
485    ! variables declaration
486    !     
487    CHARACTER(len=80),INTENT(in) :: filename
488    CHARACTER(len=1 ), OPTIONAL, INTENT(in) :: cd_type
489
490    CHARACTER*100 :: interp_type
491    CHARACTER*100 :: Child_file,Childcoordinates
492    CHARACTER*100 :: varname,childbathy
493    CHARACTER*1  :: posvar
494    CHARACTER*20,DIMENSION(:),POINTER :: Ncdf_varname
495    CHARACTER(len=20),DIMENSION(4) :: dimnames
496    !
497    REAL*8, POINTER, DIMENSION(:,:) :: lonParent,latParent => NULL()
498    REAL*8, POINTER, DIMENSION(:,:) :: lonChild,latChild,latlon_temp => NULL()
499    REAL*8, POINTER, DIMENSION(:,:,:,:) :: tabinterp4d,tabvar1,tabvar2,tabvar3 => NULL()
500    REAL*8, POINTER, DIMENSION(:,:,:) :: tabinterp3d,tabvar3d => NULL()
501    REAL*8, POINTER, DIMENSION(:) :: timedepth_temp,depth => NULL()
502    REAL*8,DIMENSION(:,:),POINTER :: matrix => NULL()
503    INTEGER,DIMENSION(:),POINTER :: src_add,dst_add  => NULL()
504    INTEGER, POINTER, DIMENSION(:) :: tabtemp1DInt => NULL()
505    REAL*8, POINTER, DIMENSION(:) :: nav_lev => NULL()
506    !     
507    LOGICAL, DIMENSION(:,:,:), POINTER :: detected_pts => NULL() 
508    REAL*8,DIMENSION(:,:,:),POINTER :: mask => NULL()
509    LOGICAL,DIMENSION(:,:),POINTER :: masksrc => NULL()
510    LOGICAL :: Interpolation,conservation,Pacifique,Extrapolation,land_level
511    !     
512    INTEGER :: deptht,time,i,status,ncid,t,ii,j,nb,numlon,numlat
513
514    !       
515    TYPE(Coordinates) :: G0,G1     
516    !
517    !*****************
518    !If coarse grid is masked possibility to activate an extrapolation process
519    !
520    Extrapolation = .FALSE.
521    PRINT*,'DEBUT INTERP_EXTRAP_VAR'
522    !
523    !*****************
524    !
525    ! check grid position
526    !
527    IF( PRESENT(cd_type) ) THEN
528       posvar = cd_type
529    ELSE
530       posvar = 'T'
531    ENDIF     
532
533    !
534    ! read dimensions in netcdf file
535    !
536    CALL Read_Ncdf_dim('x',filename,numlon)
537    CALL Read_Ncdf_dim('y',filename,numlat)
538    CALL Read_Ncdf_dim('time_counter',filename,time)
539    IF ( Dims_Existence( 'deptht' , filename ) ) THEN
540       CALL Read_Ncdf_dim('deptht',filename,deptht)
541    ELSE IF ( Dims_Existence( 'depthu' , filename ) ) THEN
542       CALL Read_Ncdf_dim('depthu',filename,deptht)
543    ELSE IF ( Dims_Existence( 'depthv' , filename ) ) THEN
544       CALL Read_Ncdf_dim('depthv',filename,deptht)
545    ELSE IF ( Dims_Existence( 'depthw' , filename ) ) THEN
546       CALL Read_Ncdf_dim('depthw',filename,deptht)
547    ELSE IF ( Dims_Existence( 'z' , filename ) ) THEN
548       CALL Read_Ncdf_dim('z',filename,deptht)
549    ELSE
550       deptht = N
551    ENDIF
552
553    !
554    ! retrieve netcdf variable name
555    !
556    CALL Read_Ncdf_VarName(filename,Ncdf_varname)
557    !
558    ! define name of child grid file
559    !
560    CALL set_child_name(filename,Child_file) 
561    CALL set_child_name(parent_coordinate_file,Childcoordinates)
562    CALL set_child_name(parent_meshmask_file,childbathy)
563    WRITE(*,*) 'Child grid file name = ',TRIM(Child_file)
564    !
565    ! create this file
566    !
567    status = nf90_create(Child_file,NF90_WRITE,ncid)
568    status = nf90_close(ncid)
569    !
570    ! read coordinates of both domains
571    !           
572    status = Read_Coordinates(parent_coordinate_file,G0)
573    status = Read_Coordinates(Childcoordinates,G1,Pacifique)
574    !
575    ! check consistency of informations read in namelist
576    !     
577    IF( imax > SIZE(G0%glamt,1) .OR. jmax > SIZE(G0%glamt,2) .OR. &
578         imax <= imin .OR. jmax <= jmin ) THEN                   
579       WRITE(*,*) 'ERROR ***** bad child grid definition ...'
580       WRITE(*,*) 'please check imin,jmin,imax,jmax,jpizoom,jpjzoom values'       
581       WRITE(*,*) ' '
582       STOP
583    ENDIF
584    !
585    IF( SIZE(G1%nav_lon,1) .NE. nxfin .OR. SIZE(G1%nav_lon,2) .NE. nyfin ) THEN
586       WRITE(*,*) 'ERROR ***** bad child coordinates or bathy file ...'
587       WRITE(*,*) ' '
588       WRITE(*,*) 'please check that child coordinates file and child bathymetry file'
589       WRITE(*,*) 'has been created with the current namelist '
590       WRITE(*,*) ' '
591       STOP
592    ENDIF
593    !     
594
595    !
596    ! Initialization of T-mask thanks to bathymetry
597    !
598    IF( Extrapolation ) THEN
599
600       WRITE(*,*) 'mask initialisation on coarse and fine grids'
601       !
602       ALLOCATE(mask(numlon,numlat,N))
603       CALL Init_mask(childbathy,G1,1,1)
604       CALL Init_mask(parent_meshmask_file,G0,1,1)
605       !
606    ENDIF
607
608    !     
609    ! select coordinates to use according to variable position
610    !
611    ALLOCATE(lonParent(numlon,numlat),latParent(numlon,numlat))
612    ALLOCATE(lonChild(nxfin,nyfin),latChild(nxfin,nyfin))
613
614    SELECT CASE(posvar)
615    CASE('T')
616       lonParent = G0%glamt
617       latParent = G0%gphit
618       lonChild  = G1%glamt
619       latChild  = G1%gphit
620       mask      = G1%tmask
621    CASE('U')
622       lonParent = G0%glamu
623       latParent = G0%gphiu
624       lonChild  = G1%glamu
625       latChild  = G1%gphiu
626       mask      = G1%umask
627    CASE('V')
628       lonParent = G0%glamv
629       latParent = G0%gphiv
630       lonChild  = G1%glamv
631       latChild  = G1%gphiv
632       mask      = G1%vmask
633    END SELECT
634
635    DEALLOCATE(G0%glamu,G0%glamv,G0%gphiu,G0%gphiv)
636    DEALLOCATE(G1%glamu,G1%glamv,G1%gphiu,G1%gphiv) 
637    DEALLOCATE(G1%glamt,G1%gphit,G0%glamt,G0%gphit)
638
639    !
640    !longitude modification if child domain covers Pacific ocean area
641    !
642    IF( lonChild(1,1) > lonChild(nxfin,nyfin) ) THEN
643       Pacifique = .TRUE.
644       WHERE( lonChild < 0 )
645          lonChild = lonChild + 360.
646       END WHERE
647       WHERE( lonParent < 0 )
648          lonParent = lonParent + 360.
649       END WHERE
650    ENDIF
651
652    !   
653    !
654    ! dimensions initialization     
655    !
656    CALL Write_Ncdf_dim('x',Child_file,nxfin)
657    CALL Write_Ncdf_dim('y',Child_file,nyfin)
658    IF ( Dims_Existence( 'deptht' , filename ) ) CALL Write_Ncdf_dim('deptht',Child_file,deptht)
659    IF ( Dims_Existence( 'depthu' , filename ) ) CALL Write_Ncdf_dim('depthu',Child_file,deptht)
660    IF ( Dims_Existence( 'depthv' , filename ) ) CALL Write_Ncdf_dim('depthv',Child_file,deptht)
661    IF ( Dims_Existence( 'depthw' , filename ) ) CALL Write_Ncdf_dim('depthw',Child_file,deptht)
662    IF ( Dims_Existence( 'z' , filename ) ) CALL Write_Ncdf_dim('z',Child_file,deptht)
663    CALL Write_Ncdf_dim('time_counter',Child_file,0)
664
665    IF( deptht .NE. 1 .AND. deptht .NE. N ) THEN
666       WRITE(*,*) '***'
667       WRITE(*,*) 'Number of vertical levels doesn t match between namelist'
668       WRITE(*,*) 'and forcing file ',TRIM(filename)
669       WRITE(*,*) 'Please check the values in namelist file'
670       WRITE(*,*) 'N = ',N
671       WRITE(*,*) 'deptht = ',deptht
672       STOP     
673    ENDIF
674    !
675    !
676    DO i = 1,SIZE(Ncdf_varname)
677       !     
678       !
679       ! ******************************LOOP ON VARIABLE NAMES*******************************************
680       !
681       !     
682       SELECT CASE (TRIM(Ncdf_varname(i)))
683          !
684          !copy nav_lon from child coordinates to output file     
685       CASE('nav_lon')   
686          ALLOCATE(latlon_temp(nxfin,nyfin))
687          CALL Read_Ncdf_var('nav_lon',TRIM(Childcoordinates),latlon_temp) 
688          CALL Write_Ncdf_var('nav_lon',(/'x','y'/),Child_file,latlon_temp,'float')
689          CALL Copy_Ncdf_att('nav_lon',TRIM(filename),Child_file, &
690               MINVAL(latlon_temp),MAXVAL(latlon_temp))
691          DEALLOCATE(latlon_temp)
692          varname = TRIM(Ncdf_varname(i))
693          Interpolation = .FALSE.
694          !     
695          !copy nav_lat from child coordinates to output file
696       CASE('nav_lat')             
697          ALLOCATE(latlon_temp(nxfin,nyfin))
698          CALL Read_Ncdf_var('nav_lat',TRIM(Childcoordinates),latlon_temp) 
699          CALL Write_Ncdf_var('nav_lat',(/'x','y'/),Child_file,latlon_temp,'float')
700          CALL Copy_Ncdf_att('nav_lat',TRIM(filename),Child_file, &
701               MINVAL(latlon_temp),MAXVAL(latlon_temp)) 
702          DEALLOCATE(latlon_temp)
703          varname = TRIM(Ncdf_varname(i))
704          Interpolation = .FALSE.
705          !
706          !copy nav_lev from restart_file to output file
707          !
708       CASE('nav_lev')
709
710          WRITE(*,*) 'copy nav_lev'
711          CALL Read_Ncdf_var('nav_lev',filename,nav_lev) 
712          IF(.NOT. dimg ) THEN
713             CALL Write_Ncdf_var('nav_lev','z',Child_file,nav_lev,'float')
714             CALL Copy_Ncdf_att('nav_lev',filename,Child_file)     
715          ENDIF
716          DEALLOCATE(nav_lev)
717          Interpolation = .FALSE.
718          !       
719          !copy time_counter from input file to output file
720       CASE('time_counter')
721          ALLOCATE(timedepth_temp(time))
722          CALL Read_Ncdf_var('time_counter',filename,timedepth_temp) 
723          CALL Write_Ncdf_var('time_counter','time_counter',  &
724               Child_file,timedepth_temp,'float')
725          CALL Copy_Ncdf_att('time_counter',TRIM(filename),Child_file)
726          DEALLOCATE(timedepth_temp)
727          varname = TRIM(Ncdf_varname(i))
728          Interpolation = .FALSE.
729          !
730          !copy deptht from input file to output file
731       CASE('deptht')
732          ALLOCATE(depth(deptht))
733          CALL Read_Ncdf_var('deptht',filename,depth) 
734          CALL Write_Ncdf_var('deptht','deptht',Child_file,depth,'float')
735          CALL Copy_Ncdf_att('deptht',TRIM(filename),Child_file)
736          varname = TRIM(Ncdf_varname(i))
737          Interpolation = .FALSE.     
738          !
739          !copy depthu from input file to output file
740       CASE('depthu')
741          ALLOCATE(depth(deptht))
742          CALL Read_Ncdf_var('depthu',filename,depth) 
743          CALL Write_Ncdf_var('depthu','depthu',Child_file,depth,'float')
744          CALL Copy_Ncdf_att('depthu',TRIM(filename),Child_file)
745          varname = TRIM(Ncdf_varname(i))
746          Interpolation = .FALSE.     
747          !
748          !copy depthv from input file to output file
749       CASE('depthv')
750          ALLOCATE(depth(deptht))
751          CALL Read_Ncdf_var('depthv',filename,depth) 
752          CALL Write_Ncdf_var('depthv','depthv',Child_file,depth,'float')
753          CALL Copy_Ncdf_att('depthv',TRIM(filename),Child_file)
754          varname = TRIM(Ncdf_varname(i))
755          Interpolation = .FALSE.     
756          !
757          !copy depthv from input file to output file
758       CASE('depthw')
759          ALLOCATE(depth(deptht))
760          CALL Read_Ncdf_var('depthw',filename,depth) 
761          CALL Write_Ncdf_var('depthw','depthw',Child_file,depth,'float')
762          CALL Copy_Ncdf_att('depthw',TRIM(filename),Child_file)
763          varname = TRIM(Ncdf_varname(i))
764          Interpolation = .FALSE.     
765          !
766          !copy time_steps from input file in case of restget use in NEMO in/out routines
767       CASE('time_steps')
768          !         print *,'avant time step'
769
770          CALL Read_Ncdf_var('time_steps',filename,tabtemp1DInt) 
771          !       print *,'timedeph = ',tabtemp1DInt
772          CALL Write_Ncdf_var('time_steps','time_counter',Child_file,tabtemp1DInt,'integer')
773          CALL Copy_Ncdf_att('time_steps',filename,Child_file) 
774          DEALLOCATE(tabtemp1DInt)
775          Interpolation = .FALSE. 
776         
777          !store tmask in output file
778       CASE('tmask')
779          dimnames(1)='x'
780          dimnames(2)='y'
781          dimnames(3)='deptht' 
782          IF (.NOT.ASSOCIATED(G1%tmask)) THEN
783             ALLOCATE(G1%tmask(nxfin,nyfin,deptht))
784             G1%tmask = 1
785          ENDIF
786          CALL Write_Ncdf_var('tmask',dimnames(1:3),Child_file,G1%tmask,'float')
787          varname = TRIM(Ncdf_varname(i))
788          Interpolation = .FALSE.     
789          !     
790       CASE default
791          varname = Ncdf_varname(i)
792          conservation = .FALSE.
793          CALL get_interptype( varname,interp_type,conservation )
794          WRITE(*,*) '**********************************************'
795          WRITE(*,*) 'varname      = ',TRIM(varname), 'at ', posvar, ' point'
796          WRITE(*,*) 'interp_type  = ',TRIM(interp_type)
797          WRITE(*,*) 'conservation = ',conservation
798          WRITE(*,*) '***********************************************'
799          Interpolation = .TRUE.
800          !         
801       END SELECT
802
803       ! //////////////// INTERPOLATION FOR 3D VARIABLES /////////////////////////////////////
804       !
805       IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 3 ) THEN
806          !     
807          ALLOCATE(detected_pts(numlon,numlat,N))
808          ALLOCATE(masksrc(numlon,numlat))                                           
809          !
810          ! ******************************LOOP ON TIME*******************************************
811          !loop on time
812          DO t = 1,time
813             !                   
814             IF(extrapolation) THEN
815                WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t
816             ELSE
817                WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t
818             ENDIF
819             !                           
820             ALLOCATE(tabvar3d(numlon,numlat,1))     
821             ALLOCATE(tabinterp3d(nxfin,nyfin,1))
822             !                   
823             CALL Read_Ncdf_var(varname,filename,tabvar3d,t)                                             
824             !
825             ! search points where extrapolation is required
826             !
827             IF(Extrapolation) THEN
828                 WHERE( tabvar3d .GE. 1.e+20 ) tabvar3d = 0.
829                IF (t .EQ. 1. ) CALL extrap_detect(G0,G1,detected_pts(:,:,1),1)
830                CALL correct_field_2d(detected_pts(:,:,1),tabvar3d,G0,masksrc,'posvar')
831             ELSE
832                masksrc = .TRUE.
833             ENDIF
834
835             IF (t.EQ.1 ) THEN
836
837                SELECT CASE(TRIM(interp_type))
838                CASE('bilinear')
839                   CALL get_remap_matrix(latParent,latChild,   &
840                        lonParent,lonChild,                    &
841                        masksrc,matrix,src_add,dst_add)
842
843                CASE('bicubic')
844                   CALL get_remap_bicub(latParent,latChild,   &
845                        lonParent,lonChild,   &
846                        masksrc,matrix,src_add,dst_add)
847
848                END SELECT
849                !                                                       
850             ENDIF
851             !     
852             SELECT CASE(TRIM(interp_type))
853             CASE('bilinear')                                                       
854                CALL make_remap(tabvar3d(:,:,1),tabinterp3d(:,:,1),nxfin,nyfin, &
855                     matrix,src_add,dst_add)     
856             CASE('bicubic')                                   
857                CALL make_bicubic_remap(tabvar3d(:,:,1),masksrc,tabinterp3d(:,:,1),nxfin,nyfin, &
858                     matrix,src_add,dst_add)                       
859             END SELECT
860             !                     
861             IF( conservation ) CALL Correctforconservation(tabvar3d(:,:,1),tabinterp3d(:,:,1), &
862                  G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 
863             !     
864             IF(Extrapolation) tabinterp3d(:,:,1) = tabinterp3d(:,:,1) * mask(:,:,1)     
865             !                     
866             dimnames(1)='x'
867             dimnames(2)='y'
868             dimnames(3)='time_counter'
869             !                                     
870             CALL Write_Ncdf_var(TRIM(varname),dimnames(1:3),&
871                  Child_file,tabinterp3d,t,'float')
872             !     
873             DEALLOCATE(tabinterp3d)
874             DEALLOCATE(tabvar3d)                     
875             !end loop on time
876          END DO
877          !                   
878          DEALLOCATE(detected_pts)
879          IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
880          DEALLOCATE( masksrc)
881
882          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)                       
883          !
884       ELSE IF( Interpolation .AND. Get_NbDims(TRIM(varname),TRIM(filename)) == 4 ) THEN
885          !
886          !
887          ! //////////////// INTERPOLATION FOR 4D VARIABLES /////////////////////////////////////
888          !
889          dimnames(1)='x'
890          dimnames(2)='y'
891          IF ( Dims_Existence( 'deptht' , filename ) ) dimnames(3)='deptht' 
892          IF ( Dims_Existence( 'depthu' , filename ) ) dimnames(3)='depthu'
893          IF ( Dims_Existence( 'depthv' , filename ) ) dimnames(3)='depthv'
894          IF ( Dims_Existence( 'depthw' , filename ) ) dimnames(3)='depthw'
895          IF ( Dims_Existence( 'z' , filename ) ) dimnames(3)='z'
896          dimnames(4)='time_counter' 
897          !
898          ! loop on vertical levels
899          !
900          DO nb = 1,deptht
901             !       
902             ALLOCATE(masksrc(numlon,numlat))
903             ALLOCATE(detected_pts(numlon,numlat,N))         
904             !
905             ! point detection et level n                                     
906             !
907             land_level = .FALSE.
908             IF( Extrapolation ) THEN
909                IF(MAXVAL(mask(:,:,nb))==0.) land_level = .TRUE.
910             ENDIF
911
912
913             IF ( land_level ) THEN
914                !
915                WRITE(*,*) 'only land points on level ',nb
916                ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
917                tabinterp4d = 0.e0
918                !                             
919                DO ii = 1,time
920                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
921                        Child_file,tabinterp4d,ii,nb,'float')
922                END DO
923                DEALLOCATE(tabinterp4d)
924                !
925             ELSE
926                !                       
927                ! loop on time                       
928                !
929                DO t = 1,time
930                   !                   
931                   ALLOCATE(tabvar1(numlon,numlat,1,1))    ! level k
932                   IF(Extrapolation) ALLOCATE(tabvar2(numlon,numlat,1,1))    ! level k-1
933                   IF(Extrapolation) ALLOCATE(tabvar3(numlon,numlat,1,1))    ! level k-2                       
934                   ALLOCATE(tabinterp4d(nxfin,nyfin,1,1)) 
935                   !
936                   IF(Extrapolation) THEN                                 
937                      IF(nb==1) THEN
938                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)   
939                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
940                      ELSE IF (nb==2) THEN
941                         CALL Read_Ncdf_var(varname,filename,tabvar2,t,nb-1)
942                         CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)           
943                         WHERE( tabvar1 .GE. 1.e+20 ) tabvar1 = 0.
944                         WHERE( tabvar2 .GE. 1.e+20 ) tabvar2 = 0.
945                      ELSE
946                         CALL Read_Ncdf_var(varname,filename,tabvar3,t,nb-2)
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                         WHERE( tabvar3 .GE. 1.e+20 ) tabvar3 = 0.
952                      ENDIF
953                      !                             
954                      IF (t.EQ.1 ) CALL extrap_detect(G0,G1,detected_pts(:,:,nb),nb)
955
956                      CALL correct_field(detected_pts(:,:,nb),tabvar1,tabvar2,&
957                           tabvar3,G0,depth,masksrc,nb,'posvar')                       
958                      DEALLOCATE(tabvar2,tabvar3)
959
960                   ELSE
961                      CALL Read_Ncdf_var(varname,filename,tabvar1,t,nb)
962                      IF(MAXVAL(tabvar1(:,:,1,1))==0.) land_level = .TRUE.
963                      masksrc = .TRUE. 
964                   ENDIF
965
966                   IF( Extrapolation ) THEN
967                      WRITE(*,*) 'interpolation/extrapolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb   
968                   ELSE
969                      WRITE(*,*) 'interpolation ',TRIM(varname),' for time t = ',t,'vertical level = ',nb
970                   ENDIF
971
972                   !                         
973
974                   IF (t.EQ.1 ) THEN
975
976                      SELECT CASE(TRIM(interp_type))
977                      CASE('bilinear')
978                         CALL get_remap_matrix(latParent,latChild,   &
979                              lonParent,lonChild,   &
980                              masksrc,matrix,src_add,dst_add)
981
982                      CASE('bicubic')
983                         CALL get_remap_bicub(latParent,latChild,   &
984                              lonParent,lonChild,   &
985                              masksrc,matrix,src_add,dst_add)
986                         !                             
987                      END SELECT
988                      !                                                       
989                   ENDIF
990                   !     
991                   SELECT CASE(TRIM(interp_type))
992                      !                             
993                   CASE('bilinear')                                                       
994                      CALL make_remap(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1),nxfin,nyfin, &
995                           matrix,src_add,dst_add)     
996                   CASE('bicubic')                                   
997                      CALL make_bicubic_remap(tabvar1(:,:,1,1),masksrc,tabinterp4d(:,:,1,1),nxfin,nyfin, &
998                           matrix,src_add,dst_add)                       
999                   END SELECT
1000                   !                     
1001                   IF( conservation ) CALL Correctforconservation(tabvar1(:,:,1,1),tabinterp4d(:,:,1,1), &
1002                        G0%e1t,G0%e2t,G1%e1t,G1%e2t,nxfin,nyfin,posvar,imin,jmin) 
1003
1004                   !
1005                   IF(Extrapolation) CALL check_extrap(G1,tabinterp4d,nb)
1006                   !     
1007                   IF(Extrapolation) tabinterp4d(:,:,1,1) = tabinterp4d(:,:,1,1) * mask(:,:,nb)     
1008                   !                     
1009                   CALL Write_Ncdf_var(TRIM(varname),dimnames, &
1010                        Child_file,tabinterp4d,t,nb,'float')
1011                   !     
1012                   DEALLOCATE(tabinterp4d)
1013                   DEALLOCATE(tabvar1)                     
1014                   !
1015                   ! end loop on time
1016                   !
1017
1018                END DO
1019
1020             ENDIF
1021
1022             !
1023             IF(ASSOCIATED(matrix)) DEALLOCATE(matrix,dst_add,src_add)
1024             DEALLOCATE( masksrc )
1025             DEALLOCATE(detected_pts)       
1026             !
1027             ! end loop on vertical levels
1028             !     
1029          END DO
1030          !   
1031          CALL Copy_Ncdf_att(TRIM(varname),TRIM(filename),Child_file)
1032          !
1033          ! fin du if interpolation ...
1034          !
1035       ENDIF
1036       !
1037    END DO
1038
1039    PRINT *,'FIN DE INTERPEXTRAPVAR'
1040    !   
1041    IF(Extrapolation) DEALLOCATE(G1%tmask,G0%tmask)
1042    DEALLOCATE(G0%e1t,G0%e2t,G1%e1t,G1%e2t)
1043    IF(ASSOCIATED(depth)) DEALLOCATE(depth)   
1044    !
1045  END SUBROUTINE Interp_Extrap_var
1046  !
1047  !
1048END MODULE agrif_interpolation
Note: See TracBrowser for help on using the repository browser.