source: utils/tools/NESTING/src/agrif_extrapolation.f90 @ 10383

Last change on this file since 10383 was 10383, checked in by clem, 23 months ago

ice restart should work in the nesting tools now. However ocean restart has been broken for some time

  • Property svn:keywords set to Id
File size: 20.6 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_extrapolation
9  !
10  USE agrif_types
11  USE agrif_readwrite 
12  USE io_netcdf 
13  USE agrif_gridsearch 
14 
15  IMPLICIT NONE
16
17CONTAINS
18  !
19  !************************************************************************
20  !                           *
21  ! MODULE  AGRIF_EXTRAPOLATION                 *
22  !                           *
23  !************************************************************************   
24  !     
25  !****************************************************************
26  !     subroutine extrap_detect             *
27  !                        *
28  !     detection on each level of points          *
29  !     where extrapolation is required            *
30  !                        *
31  !****************************************************************         
32  !     
33  !
34  SUBROUTINE extrap_detect(G0,G1,detected,n,posvar) 
35    !     
36    LOGICAL, DIMENSION(:,:) :: detected
37    TYPE(Coordinates) :: G0,G1       
38    CHARACTER(*), OPTIONAL :: posvar
39    INTEGER :: i,j,k,ic,jc,compt,dst_add,n
40    INTEGER, DIMENSION(1) :: i_min,j_min     
41    !                               
42    IF( PRESENT(posvar) .AND. posvar == 'U' ) THEN     
43       CALL get_detected_pts(G0%gphiu,G1%gphiu,G0%glamu,G1%glamu,   &
44            G0%umask(:,:,n),G1%umask(:,:,n),detected(:,:))       
45    ELSEIF( PRESENT(posvar) .AND. posvar == 'V' ) THEN
46       CALL get_detected_pts(G0%gphiv,G1%gphiv,G0%glamv,G1%glamv,   &
47            G0%vmask(:,:,n),G1%vmask(:,:,n),detected(:,:))                                 
48    ELSEIF( PRESENT(posvar) .AND. posvar == 'F' ) THEN
49       CALL get_detected_pts(G0%gphif,G1%gphif,G0%glamf,G1%glamf,   &
50            G0%fmask(:,:,n),G1%fmask(:,:,n),detected(:,:))                                 
51    ELSE
52       CALL get_detected_pts(G0%gphit,G1%gphit,G0%glamt,G1%glamt,   &
53            G0%tmask(:,:,n),G1%tmask(:,:,n),detected(:,:))                                 
54!clem       CALL get_detected_pts(G0%nav_lat,G1%nav_lat,G0%nav_lon,G1%nav_lon,   &
55!clem            G0%tmask(:,:,n),G1%tmask(:,:,n),detected(:,:))       
56    ENDIF
57    !     
58  END SUBROUTINE extrap_detect
59  !     
60  !   
61  !****************************************************************
62  !     end subroutine extrap_detect            *
63  !****************************************************************
64  !     
65  !     
66  !****************************************************************
67  !    subroutine correct_field              *
68  ! correct field on detected points            *
69  !                        *
70  !****************************************************************         
71  !
72  SUBROUTINE correct_field(detected_pts,tabin,tabinkm1,tabinkm2,G0,nav_lev,newmask,k,posvar)
73    !
74    LOGICAL, DIMENSION(:,:) :: detected_pts
75    LOGICAL, DIMENSION(:,:) :: newmask
76    CHARACTER(*),OPTIONAL :: posvar
77    INTEGER :: k
78    !
79    INTEGER :: i,j,ii,jj,nx,ny,n,lbi,ubi,lbj,ubj,kpos,ipos,jpos,r
80    !
81    REAL*8, DIMENSION(:,:,:,:) ::  tabin
82    REAL*8, DIMENSION(:,:,:,:) ::  tabinkm1
83    REAL*8, DIMENSION(:,:,:,:) ::  tabinkm2
84    REAL*8, DIMENSION(:) :: nav_lev
85    REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: mask
86    REAL*8, DIMENSION(:,:), ALLOCATABLE :: lon,lat 
87    REAL*8 :: deriv,deriv_min
88    LOGICAL :: found
89    !     
90    TYPE(Coordinates) :: G0 
91    !
92    ! copy coarse grid mask in newmask
93    !           
94    IF ( PRESENT(posvar) .AND. posvar == 'U' ) THEN
95       WHERE(G0%umask(:,:,k) == 1. )
96          newmask(:,:) = .TRUE.
97       ELSEWHERE
98          newmask(:,:) = .FALSE.
99       END WHERE
100       ALLOCATE(mask(SIZE(G0%umask,1),SIZE(G0%umask,2),SIZE(G0%umask,3)))
101       ALLOCATE(lat(SIZE(G0%umask,1),SIZE(G0%umask,2)))
102       ALLOCATE(lon(SIZE(G0%umask,1),SIZE(G0%umask,2)))
103       mask = G0%umask
104       lat = G0%gphiu 
105       lon = G0%glamu     
106    ELSE IF ( PRESENT(posvar) .AND. posvar == 'V' ) THEN
107       WHERE(G0%vmask(:,:,k) == 1. )
108          newmask(:,:) = .TRUE.
109       ELSEWHERE
110          newmask(:,:) = .FALSE.
111       END WHERE
112       ALLOCATE(mask(SIZE(G0%vmask,1),SIZE(G0%vmask,2),SIZE(G0%vmask,3)))
113       ALLOCATE(lat(SIZE(G0%vmask,1),SIZE(G0%vmask,2)))
114       ALLOCATE(lon(SIZE(G0%vmask,1),SIZE(G0%vmask,2)))
115       mask = G0%vmask       
116       lat = G0%gphiv 
117       lon = G0%glamv
118    ELSE IF ( PRESENT(posvar) .AND. posvar == 'F' ) THEN
119       WHERE(G0%fmask(:,:,k) == 1. )
120          newmask(:,:) = .TRUE.
121       ELSEWHERE
122          newmask(:,:) = .FALSE.
123       END WHERE
124       ALLOCATE(mask(SIZE(G0%fmask,1),SIZE(G0%fmask,2),SIZE(G0%fmask,3)))
125       ALLOCATE(lat(SIZE(G0%fmask,1),SIZE(G0%fmask,2)))
126       ALLOCATE(lon(SIZE(G0%fmask,1),SIZE(G0%fmask,2)))
127       mask = G0%fmask       
128       lat = G0%gphif 
129       lon = G0%glamf
130    ELSE     
131       WHERE(G0%tmask(:,:,k) == 1. )
132          newmask(:,:) = .TRUE.
133       ELSEWHERE
134          newmask(:,:) = .FALSE.
135       END WHERE
136       ALLOCATE(mask(SIZE(G0%tmask,1),SIZE(G0%tmask,2),SIZE(G0%tmask,3)))
137       ALLOCATE(lat(SIZE(G0%tmask,1),SIZE(G0%tmask,2)))
138       ALLOCATE(lon(SIZE(G0%tmask,1),SIZE(G0%tmask,2)))
139       mask = G0%tmask       
140       lat = G0%gphit 
141       lon = G0%glamt
142!clem       mask = G0%tmask
143!clem       lon = G0%nav_lon
144!clem       lat = G0%nav_lat
145    ENDIF
146    !
147    ! dimensions initialisation
148    !
149    nx = SIZE(tabin,1)
150    ny = SIZE(tabin,2)   
151    !
152    !       
153    DO j = 1,ny   
154       !
155       DO i = 1,nx         
156          !
157          IF( detected_pts(i,j) ) THEN       
158             !
159             r = 0
160             found = .FALSE.
161             deriv_min = 2000000.
162             ipos=0
163             jpos=0
164             !
165             DO WHILE(.NOT. found)
166                !
167                r = r + 1     
168                !       
169                IF(i-r < 1 ) THEN
170                   lbi = 1
171                   ubi = MIN(i+r,nx)
172                ELSE IF(i+r > nx) THEN
173                   lbi = MAX(i-r,1)
174                   ubi = nx
175                ELSE
176                   lbi = i-r
177                   ubi = i+r
178                ENDIF
179                !
180                IF(j-r < 1) THEN
181                   lbj = 1
182                   ubj = MIN(j+r,ny)
183                ELSE IF(j+r > ny) THEN
184                   lbj = MAX(j-r,1)
185                   ubj = ny
186                ELSE
187                   lbj = j-r
188                   ubj = j+r
189                ENDIF
190                !
191                DO jj = lbj,ubj,ubj-lbj
192                   DO ii = lbi,ubi,ubi-lbi
193                      !
194                      deriv = search_pts_h(ii,jj,k,i,j,k,tabin(:,:,1,1),mask,lon,lat) 
195                      !
196                      IF( ABS(deriv) < deriv_min ) THEN                                                   
197                         deriv_min = ABS(deriv)
198                         ipos = ii
199                         jpos = jj
200                         kpos = k
201                      ENDIF
202                      !                           
203                      deriv = search_pts_v(ii,jj,k-1,i,j,k,tabinkm1,tabinkm2,mask,nav_lev,lon,lat)
204                      !
205                      IF( ABS(deriv) < deriv_min ) THEN                                                   
206                         deriv_min = ABS(deriv)
207                         ipos = ii
208                         jpos = jj
209                         kpos = k-1
210                      ENDIF
211                      !                                       
212                   END DO
213                END DO
214                !
215                !                                 
216                IF( deriv_min < 2000000.  ) THEN 
217                   !               
218                   IF(kpos == k)   tabin(i,j,1,1) = tabin(ipos,jpos,1,1)
219                   IF(kpos == k-1) tabin(i,j,1,1) = tabinkm1(ipos,jpos,1,1)   
220                   found = .TRUE.
221                   newmask(i,j) = .TRUE.
222                ELSE IF ((lbi == 1).AND.(ubi == nx).AND.(lbj == 1).AND.(ubj == ny)) THEN
223                   found = .TRUE.
224                   newmask(i,j) = .FALSE.
225                   !
226                ENDIF
227                !
228             END DO !do while
229             !
230          ENDIF
231          !
232       END DO
233       !
234    END DO
235    !
236    DEALLOCATE(mask,lon,lat)
237    !
238  END SUBROUTINE correct_field
239  !     
240  !**************************************************************
241  !    end subroutine correct_field
242  !**************************************************************         
243 
244  SUBROUTINE correct_field_2d(detected_pts,tabin,G0,newmask,posvar)
245    !
246    LOGICAL, DIMENSION(:,:) :: detected_pts
247    LOGICAL, DIMENSION(:,:) :: newmask
248    CHARACTER(*), OPTIONAL :: posvar
249    LOGICAL :: found
250    INTEGER :: k
251    !
252    INTEGER :: i,j,ii,jj,nx,ny,n,lbi,ubi,lbj,ubj,ipos,jpos,r
253    !
254    REAL*8, DIMENSION(:,:,:), ALLOCATABLE :: mask
255    REAL*8, DIMENSION(:,:), ALLOCATABLE :: lon,lat 
256    REAL*8, DIMENSION(:,:,:) ::  tabin
257    REAL*8 :: deriv,deriv_min
258    !     
259    TYPE(Coordinates) :: G0 
260    !
261    ! copy coarse grid mask in newmask
262    !           
263    mask = G0%tmask
264    lon = G0%nav_lon
265    lat = G0%nav_lat     
266    IF ( PRESENT(posvar) .AND. posvar == 'U' ) THEN
267       WHERE(G0%umask(:,:,1) == 1. )
268          newmask(:,:) = .TRUE.
269       ELSEWHERE
270          newmask(:,:) = .FALSE.
271       END WHERE
272       ALLOCATE(mask(SIZE(G0%umask,1),SIZE(G0%umask,2),SIZE(G0%umask,3)))
273       ALLOCATE(lat(SIZE(G0%umask,1),SIZE(G0%umask,2)))
274       ALLOCATE(lon(SIZE(G0%umask,1),SIZE(G0%umask,2)))
275       mask = G0%umask
276       lat = G0%gphiu 
277       lon = G0%glamu     
278    ELSE IF ( PRESENT(posvar) .AND. posvar == 'V' ) THEN
279       WHERE(G0%vmask(:,:,1) == 1. )
280          newmask(:,:) = .TRUE.
281       ELSEWHERE
282          newmask(:,:) = .FALSE.
283       END WHERE
284       ALLOCATE(mask(SIZE(G0%vmask,1),SIZE(G0%vmask,2),SIZE(G0%vmask,3)))
285       ALLOCATE(lat(SIZE(G0%vmask,1),SIZE(G0%vmask,2)))
286       ALLOCATE(lon(SIZE(G0%vmask,1),SIZE(G0%vmask,2)))
287       mask = G0%vmask       
288       lat = G0%gphiv 
289       lon = G0%glamv
290    ELSE     
291       WHERE(G0%tmask(:,:,1) == 1. )
292          newmask(:,:) = .TRUE.
293       ELSEWHERE
294          newmask(:,:) = .FALSE.
295       END WHERE
296       ALLOCATE(mask(SIZE(G0%tmask,1),SIZE(G0%tmask,2),SIZE(G0%tmask,3)))
297       ALLOCATE(lat(SIZE(G0%tmask,1),SIZE(G0%tmask,2)))
298       ALLOCATE(lon(SIZE(G0%tmask,1),SIZE(G0%tmask,2)))
299       mask = G0%tmask
300       lon = G0%nav_lon
301       lat = G0%nav_lat
302    ENDIF
303
304    !
305    ! dimensions initialisation
306    !
307    nx = SIZE(tabin,1)
308    ny = SIZE(tabin,2)   
309    !       
310    DO i = 1,nx         
311       !
312       DO j = 1,ny   
313          !                   
314          !                     
315          IF( detected_pts(i,j) ) THEN       
316             !
317             r = 0 
318             found = .FALSE.
319             deriv_min = 2000000.
320             ipos=0
321             jpos=0
322             !
323             DO WHILE (.NOT. found )
324
325                !
326                r = r + 1 
327                !
328                IF(i-r < 1 ) THEN
329                   lbi = 1
330                   ubi = MIN(i+r,nx)
331                ELSE IF(i+r > nx) THEN
332                   lbi = MAX(i-r,1)
333                   ubi = nx
334                ELSE
335                   lbi = i-r
336                   ubi = i+r
337                ENDIF
338                !
339                IF(j-r < 1) THEN
340                   lbj = 1
341                   ubj = MIN(j+r,ny)
342                ELSE IF(j+r > ny) THEN
343                   lbj = MAX(j-r,1)
344                   ubj = ny
345                ELSE
346                   lbj = j-r
347                   ubj = j+r
348                ENDIF
349                !                                 
350                DO ii = lbi,ubi
351                   DO jj = lbj,ubj
352                      !
353                      deriv = search_pts_h(ii,jj,1,i,j,1,tabin(:,:,1),mask,lon,lat) 
354                      !
355                      IF( ABS(deriv) < deriv_min ) THEN                                                   
356                         deriv_min = ABS(deriv)
357                         ipos = ii
358                         jpos = jj
359                      ENDIF
360                      !                                       
361                   END DO
362                END DO
363                !
364                !                                   
365                IF( deriv_min < 2000000.  ) THEN 
366                   !
367                   found = .TRUE.                                             
368                   tabin(i,j,1) = tabin(ipos,jpos,1) 
369                   newmask(i,j) = .TRUE.
370                   !
371                ENDIF
372                !
373             END DO !do while
374             !
375          ENDIF
376          !
377       END DO
378       !
379    END DO
380    !
381    DEALLOCATE(mask,lon,lat)
382    !
383  END SUBROUTINE correct_field_2d
384  !     
385  !**************************************************************
386  !    function get_dist
387  !**************************************************************         
388  !
389
390  !
391  REAL*8 FUNCTION  get_dist(plat1,plon1,plat2,plon2)
392    !
393    REAL*8 :: plat1,plon1,plat2,plon2
394    REAL*8 :: dist,ra,rad,rpi,lat,lon
395    !     
396    rpi = 3.141592653589793
397    rad = rpi/180.
398    ra  = 6371229.   
399    !     
400    lat = plat2-plat1
401    lon = plon2-plon1
402    !
403    dist = ra * rad * SQRT( (COS(rad*(plat1+plat2)/2.)*lon)**2 + lat**2 )         
404    get_dist = dist
405    RETURN
406    !
407  END FUNCTION get_dist
408
409  !
410  !     
411  !**************************************************************
412  !    end function get_dist
413  !**************************************************************
414  !
415  !     
416  !**************************************************************
417  !    subroutine check_extrap
418  !**************************************************************         
419  !
420
421  !
422  SUBROUTINE check_extrap(Grid,tabin,k)
423    !
424    REAL*8, DIMENSION(:,:,:,:) ::  tabin
425    TYPE(Coordinates) :: Grid 
426    INTEGER :: i,j,k 
427   
428    DO i = 2,SIZE(tabin,1)-1
429       DO j=2,SIZE(tabin,2)-1
430          !                     
431          IF( Grid%tmask(i,j,k) == 1. .AND. tabin(i,j,1,1)==0.) THEN
432             !     
433             WRITE(*,*) 'no masked point with value zero (',i,',',j,',',k,')'
434             !                   
435          ENDIF
436
437       END DO
438    END DO
439    !
440  END SUBROUTINE check_extrap
441
442  !
443  !     
444  !**************************************************************
445  !    end subroutine check_extrap
446  !************************************************************** 
447  !
448  !**************************************************************
449  !    subroutine search_pts_h
450  !**************************************************************
451  !
452  REAL*8 FUNCTION search_pts_h(i,j,k,ipt,jpt,kpt,tabvar,mask,lon,lat)
453    !
454    REAL*8 :: hx,hy,fx,fy
455    REAL*8 :: h_x,h_y
456    REAL*8, DIMENSION(:,:) :: tabvar
457    INTEGER :: i,j,k,ipt,jpt,kpt,nx,ny
458    LOGICAL :: foundx,foundy
459    REAL*8, DIMENSION(:,:,:) :: mask
460    REAL*8, DIMENSION(:,:) :: lon,lat
461    !
462    !
463    foundx = .TRUE.
464    foundy = .TRUE.
465    !     
466    nx = SIZE(tabvar,1)
467    ny = SIZE(tabvar,2)
468    !
469    IF( i==ipt .AND. j==jpt ) THEN   
470       search_pts_h = 2000000.
471       RETURN
472    ENDIF
473    !     
474    IF( mask(i,j,k) == 0. ) THEN   
475       search_pts_h = 2000000.
476       RETURN
477    ENDIF
478    !
479    ! x direction
480    !
481    IF(i+1<=nx .AND. i-1>=1) THEN
482       IF(mask(i+1,j,k)==1. .AND. mask(i-1,j,k)==1.) THEN     
483          hx = get_dist(lat(i+1,j),lon(i+1,j),&
484               lat(i-1,j),lon(i-1,j))
485          fx = (tabvar(i+1,j) - tabvar(i-1,j))/hx
486       ELSE IF(mask(i+1,j,k)==1. .AND. mask(i-1,j,k)==0. .AND. mask(i,j,k)==1.) THEN
487          hx = get_dist(lat(i+1,j),lon(i+1,j),&
488               lat(i,j),lon(i,j))
489          fx = (tabvar(i+1,j) - tabvar(i,j))/hx
490       ELSE IF(mask(i+1,j,k)==0. .AND. mask(i-1,j,k)==1. .AND. mask(i,j,k)==1.) THEN
491          hx = get_dist(lat(i,j),lon(i,j),&
492               lat(i-1,j),lon(i-1,j))
493          fx = (tabvar(i,j) - tabvar(i-1,j))/hx
494       ELSE
495          foundx = .FALSE.                 
496       ENDIF
497       !           
498    ELSE IF(i+1<=nx .AND. i>=1) THEN   
499       !
500       IF(mask(i+1,j,k)==1. .AND. mask(i,j,k)==1.) THEN     
501          hx = get_dist(lat(i+1,j),lon(i+1,j),&
502               lat(i,j),lon(i,j))
503          fx = (tabvar(i+1,j) - tabvar(i,j))/hx
504       ELSE
505          foundx = .FALSE.             
506       ENDIF
507       !   
508    ELSE IF(i<=nx .AND. i-1>=1) THEN   
509       !
510       IF(mask(i,j,k)==1. .AND. mask(i-1,j,k)==1.) THEN     
511          hx = get_dist(lat(i,j),lon(i,j),&
512               lat(i-1,j),lon(i-1,j))
513          fx = (tabvar(i,j) - tabvar(i-1,j))/hx
514       ELSE
515          foundx = .FALSE.           
516       ENDIF
517       !
518    ELSE
519        foundy = .FALSE.             
520    ENDIF
521
522    !
523    ! y direction
524    !
525    IF(j+1<=ny .AND. j-1>=1) THEN     
526       IF( mask(i,j+1,k)==1. .AND. mask(i,j-1,k)==1. ) THEN     
527          hy = get_dist(lat(i,j+1),lon(i,j+1),&
528               lat(i,j-1),lon(i,j-1))                     
529          fy = (tabvar(i,j+1) - tabvar(i,j-1))/hy
530       ELSE IF( mask(i,j+1,k)==1. .AND. mask(i,j-1,k)==0. .AND. mask(i,j,k)==1.) THEN     
531          hy = get_dist(lat(i,j+1),lon(i,j+1),&
532               lat(i,j),lon(i,j))                     
533          fy = (tabvar(i,j+1) - tabvar(i,j))/hy     
534       ELSE IF( mask(i,j+1,k)==0. .AND. mask(i,j-1,k)==1. .AND. mask(i,j,k)==1.) THEN     
535          hy = get_dist(lat(i,j),lon(i,j),&
536               lat(i,j-1),lon(i,j-1))                     
537          fy = (tabvar(i,j) - tabvar(i,j-1))/hy     
538       ELSE
539          foundy = .FALSE.                   
540       ENDIF
541       !           
542    ELSE IF(j+1<=ny .AND. j>=1) THEN   
543       !
544       IF(mask(i,j+1,k)==1. .AND. mask(i,j,k)==1.) THEN     
545          hy = get_dist(lat(i,j+1),lon(i,j+1),&
546               lat(i,j),lon(i,j))
547          fy = (tabvar(i,j+1) - tabvar(i,j))/hy
548       ELSE
549          foundy = .FALSE.           
550       ENDIF
551       !   
552    ELSE IF(j<=ny .AND. j-1>=1) THEN   
553       !
554       IF(mask(i,j,k)==1. .AND. mask(i,j-1,k)==1.) THEN     
555          hy = get_dist(lat(i,j),lon(i,j),&
556               lat(i,j-1),lon(i,j-1))
557          fy = (tabvar(i,j) - tabvar(i,j-1))/hy
558       ELSE
559          foundy = .FALSE.             
560       ENDIF
561       !
562    ELSE
563        foundy = .FALSE.             
564    ENDIF
565    !                   
566    h_x = get_dist(lat(ipt,jpt),lon(ipt,jpt),lat(ipt,jpt),lon(i,j))
567    h_y = get_dist(lat(ipt,jpt),lon(ipt,jpt),lat(i,j),lon(ipt,jpt))
568    !
569    IF(.NOT.foundx .AND. .NOT.foundy)THEN     
570       search_pts_h = 2000000.
571    ELSE IF( foundx .AND. foundy) THEN       
572       search_pts_h = h_x * fx + h_y * fy
573    ELSE IF( .NOT.foundx .AND. foundy .AND. h_y.NE.0.) THEN       
574       search_pts_h = h_y * fy
575    ELSE IF( foundx .AND. .NOT.foundy .AND. h_x.NE.0.) THEN       
576       search_pts_h = h_x * fx
577    ELSE   
578       search_pts_h = 2000000.             
579    ENDIF
580
581    !     
582    RETURN         
583    !
584  END FUNCTION search_pts_h
585  !
586  !**************************************************************
587  !    end subroutine search_pts_h
588  !**************************************************************
589  !
590  !**************************************************************
591  !    subroutine search_pts_v
592  !**************************************************************
593  !
594  REAL*8 FUNCTION search_pts_v(i,j,k,ipt,jpt,kpt,tabvarkm1,tabvarkm2,mask,depth,lon,lat)
595    !
596    REAL*8 :: hz,fz,dz,fh
597    REAL*8, DIMENSION(:) :: depth 
598    REAL*8, DIMENSION(:,:,:,:) :: tabvarkm1,tabvarkm2
599    INTEGER :: i,j,k,ipt,jpt,kpt,nx,ny
600    LOGICAL :: foundz
601    REAL*8, DIMENSION(:,:,:) :: mask
602    REAL*8, DIMENSION(:,:) :: lon,lat
603    !         
604    IF( k <= 2 .OR. mask(i,j,k) == 0.   ) THEN
605       !
606       search_pts_v = 2000000.
607       RETURN
608       !
609    ELSE IF( i==ipt .AND. j==jpt .AND. mask(i,j,k-1) == 1. .AND. mask(i,j,k-2) == 1. ) THEN
610       !
611       dz = depth(k) - depth(k-1)
612       hz = depth(kpt) - depth(k)     
613       search_pts_v = ((tabvarkm2(i,j,1,1) - tabvarkm1(i,j,1,1))/dz)*hz
614       RETURN
615       !
616    ELSE
617       !
618       IF( mask(i,j,k) == 1. .AND. mask(i,j,k-1) == 1. ) THEN
619          !                         
620          dz = depth(k) - depth(k-1)
621          fz = (tabvarkm2(i,j,1,1) - tabvarkm1(i,j,1,1))/dz
622          hz = depth(kpt) - depth(k)
623          !
624       ELSE
625          foundz = .FALSE.
626       ENDIF
627       !             
628       fh = search_pts_h(i,j,k,ipt,jpt,k,tabvarkm1(:,:,1,1),mask,lon,lat)
629       !
630       IF(foundz) THEN
631          search_pts_v = hz * fz + fh
632          RETURN
633       ELSE       
634          search_pts_v = 2000000. 
635          RETURN
636       ENDIF
637       !
638    ENDIF
639    WRITE(*,*) 'cas 2', search_pts_v
640    !     
641    RETURN         
642    !
643  END FUNCTION search_pts_v
644  !
645  !**************************************************************
646  !    end subroutine search_pts_v
647  !**************************************************************
648  !
649  !
650END MODULE agrif_extrapolation
Note: See TracBrowser for help on using the repository browser.