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

Last change on this file since 12253 was 10383, checked in by clem, 2 years 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.