source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/TOOLS/NESTING/src/agrif_interpolation.f90 @ 5445

Last change on this file since 5445 was 5445, checked in by davestorkey, 5 years ago

Clear SVN keywords from 2015/dev_r5021_UKMO1_CICE_coupling branch.

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