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

Last change on this file since 12253 was 12253, checked in by clem, 10 months ago

1) resolve some conflicts when using partial steps or not. 2) Make sure grids remain identical when they should. 3) in case of a double zoom, it is no more necessary to have the bilinear interpolation option when updating the 1st parent grid. Btw, I know these comments are unclear but it is very difficult to explain…

  • Property svn:keywords set to Id
File size: 28.7 KB
Line 
1!************************************************************************
2! Fortran 95 OPA Nesting tools                  *
3!                          *
4!     Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr)   *
5!                        Laurent Debreu (Laurent.Debreu@imag.fr)  *
6!************************************************************************
7!
8!Smoothing procedures : Pierrick Penven 2004
9!
10MODULE agrif_connect_topo
11  !
12  USE agrif_types
13  !
14  IMPLICIT NONE           
15  !
16CONTAINS
17  !
18  !
19  !************************************************************************
20  !                           *
21  ! MODULE  CONNECT_TOPO                     *
22  !                           *
23  ! module containing subroutine used for :           *
24  !   - Parent-Child bathymetry connection            *
25  !   - Bathymetry smoothing                 *
26  !   - Meters to levels conversion             *
27  !   - Parent Bathymetry update                   *
28  !                           *
29  !************************************************************************
30  !
31  !****************************************************************
32  !   subroutine init_constant_bathy            *
33  !                        *
34  !                        *
35  ! - input :                    *
36  !     coarse_bathy : coarse grid bathymetry         *
37  ! - ouput :                    *
38  !    bathy_fin_constant : coarse bathymetry on fine grid  *
39  !                        *
40  !****************************************************************
41  !
42  SUBROUTINE init_constant_bathy(coarse_bathy,bathy_fin_constant)
43    !       
44    IMPLICIT NONE 
45    !     
46    INTEGER :: i,j,ii,jj,ji
47    INTEGER :: jpt,ipt,diff,indx,indy,bornex,borney,bornex2,borney2
48    INTEGER :: jdeb,ideb,ifin,jfin
49    REAL*8, DIMENSION(:,:)  :: coarse_bathy
50    REAL*8, DIMENSION(:,:),POINTER :: bathy_fin_constant
51    TYPE(Coordinates) :: Grid     
52    !
53    IF( ln_agrif_domain ) THEN
54    !
55       diff = 0
56       IF(MOD(rho,2) .EQ. 0) diff = 1
57       !
58       indx = 1 + nbghostcellsfine + CEILING(irafx/2.0) + diff     
59       indy = 1 + nbghostcellsfine + CEILING(irafy/2.0) + diff
60       bornex = 1+nbghostcellsfine + CEILING(irafx/2.0) + diff - irafx
61       borney = 1+nbghostcellsfine + CEILING(irafy/2.0) + diff - irafy
62       bornex2 = nxfin - (nbghostcellsfine) + irafx - CEILING(irafx/2.0) 
63       borney2 = nyfin - (nbghostcellsfine) + irafy - CEILING(irafy/2.0) 
64       !   
65       ALLOCATE(bathy_fin_constant(bornex-FLOOR(irafx/2.0):bornex2+FLOOR(irafx/2.0), &
66          borney-FLOOR(irafy/2.0):borney2+FLOOR(irafy/2.0)))
67       !
68       DO j = borney,borney2,irafy
69         
70          jpt = jmin + 1 + nbghostcellsfine + (j-indy)/irafy
71          IF(j<=1) jpt = jmin + 1
72         
73          DO i = bornex,bornex2,irafx
74             
75             ipt = imin + 1 + nbghostcellsfine + (i-indx)/irafx
76             IF(i<=1) ipt = imin + 1
77             !       
78             DO jj = j-FLOOR(irafy/2.0),j+FLOOR(irafy/2.0)-diff
79                DO ii = i-FLOOR(irafx/2.0),i+FLOOR(irafx/2.0)-diff
80                   
81                   bathy_fin_constant(ii,jj) = coarse_bathy(ipt,jpt)
82                   
83                END DO
84             END DO
85             
86          END DO
87       END DO
88       
89    ELSE
90
91       ALLOCATE(bathy_fin_constant(1:nxfin,1:nyfin))
92
93       DO j = 1,nyfin-irafy+1,irafy
94          jpt = jmin + FLOOR( (j - 1.) / irafy )
95          !
96          DO i = 1,nxfin-irafx+1,irafx
97             ipt = imin + FLOOR( (i - 1.) / irafx )
98             !
99             bathy_fin_constant(i:i+irafx-1,j:j+irafy-1) = coarse_bathy(ipt,jpt)
100             !
101          END DO
102       END DO
103 
104    ENDIF
105    !
106    !
107  END SUBROUTINE init_constant_bathy
108  !
109  !****************************************************************
110  !   subroutine meter_to_levels             *
111  !                        *
112  ! subroutine to convert bathymetry in meters to bathymetry   *
113  ! in vertical levels                 *
114  !                        *
115  ! - input/output :                *
116  !     Grid : grid where conversion is required         *
117  !                        *
118  !various input parameters come from namelist.input files  *
119  !****************************************************************
120  !
121  SUBROUTINE meter_to_levels(Grid)
122    !
123    IMPLICIT NONE
124    !
125    REAL*8 :: za2,za1,za0,zsur,zacr,zkth,zacr2,zkth2,zmin,zmax
126    TYPE(Coordinates) :: Grid
127    INTEGER :: i,j
128    INTEGER, DIMENSION(1) :: k
129    INTEGER :: k1,ji,jj,jpi,jpj
130    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
131    !       
132    WRITE(*,*) 'convert bathymetry from etopo to vertical levels'
133    !
134    jpi = SIZE(Grid%bathy_meter,1)
135    jpj = SIZE(Grid%bathy_meter,2)
136    !                     
137    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
138         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
139       !   
140       za1=( ppdzmin - pphmax / (N-1) )          &
141            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
142            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
143            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
144
145       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
146       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
147       !
148    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
149         pa0.NE.0 .AND. pa1.NE.0 ) THEN
150       !       
151       zsur = psur
152       za0  = pa0
153       za1  = pa1
154       za2  = pa2
155       !
156    ELSE
157       !       
158       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
159       WRITE(*,*) ' '
160       WRITE(*,*) 'please check values of variables'
161       WRITE(*,*) 'in namelist vertical_grid section'
162       WRITE(*,*) ' ' 
163       STOP     
164       !       
165    ENDIF
166
167    zacr = ppacr
168    zkth = ppkth
169    zacr2 = ppacr2
170    zkth2 = ppkth2   
171    !
172    ALLOCATE(gdepw(N),gdept(N),e3w(N),e3t(N))
173    !
174    IF( ppkth == 0. ) THEN            !  uniform vertical grid
175       za1 = pphmax / FLOAT(N-1) 
176       DO i = 1, N
177          gdepw(i) = ( i - 1   ) * za1
178          gdept(i) = ( i - 0.5 ) * za1
179          e3w  (i) =  za1
180          e3t  (i) =  za1
181       END DO
182    ELSE                            ! Madec & Imbard 1996 function
183       IF( .NOT. ldbletanh ) THEN
184          DO i = 1,N
185             !
186             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
187             gdept(i) = (zsur+za0*(i+0.5)+za1*zacr*LOG(COSH(((i+0.5)-zkth)/zacr)))
188             e3w(i)   = (za0 + za1 * TANH((i-zkth)/zacr))
189             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
190             !
191          END DO
192       ELSE
193          DO i = 1,N
194             ! Double tanh function
195             gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
196                &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
197             gdept(i) = ( zsur + za0*(i+0.5) + za1 * zacr * LOG ( COSH( ((i+0.5)-zkth ) / zacr  ) )    &
198                &                            + za2 * zacr2* LOG ( COSH( ((i+0.5)-zkth2) / zacr2 ) )  )
199             e3w  (i) =          za0         + za1        * TANH(       (i-zkth ) / zacr  )            &
200                &                            + za2        * TANH(       (i-zkth2) / zacr2 )
201             e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
202                &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
203            END DO
204       ENDIF
205    ENDIF
206    !
207    gdepw(1) = 0.0
208    !
209    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
210       !
211       DO i = 1, N-1 
212          e3t(i) = gdepw(i+1)-gdepw(i)
213       END DO
214       e3t(N) = e3t(N-1)
215
216       DO i = 2, N 
217          e3w(i) = gdept(i) - gdept(i-1)
218       END DO
219       e3w(1  ) = 2. * (gdept(1) - gdepw(1))
220    END IF 
221    !
222    zmax = gdepw(N) + e3t(N)
223    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
224    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
225    ENDIF
226    zmin = gdepw(i+1)
227    !
228    IF ( .NOT. ASSOCIATED(Grid%bathy_level)) &
229         ALLOCATE(Grid%bathy_level(jpi,jpj))   
230    !
231    Grid%bathy_level = N-1
232    !
233    DO jj = 1, jpj
234       DO ji= 1, jpi
235          IF( Grid%bathy_meter(ji,jj) <= 0. )   &
236               Grid%bathy_level(ji,jj) = INT( Grid%bathy_meter(ji,jj) )
237       END DO
238    END DO
239    !
240    DO jj = 1, jpj
241       DO ji= 1, jpi
242          IF( Grid%bathy_meter(ji,jj) <= 0. ) THEN
243             Grid%bathy_meter(ji,jj) = 0.e0
244          ELSE
245             Grid%bathy_meter(ji,jj) = MAX( Grid%bathy_meter(ji,jj), zmin )
246             Grid%bathy_meter(ji,jj) = MIN( Grid%bathy_meter(ji,jj), zmax )
247          ENDIF
248       END DO
249    END DO
250    !
251    !
252    !
253    DO j = 1,jpj
254       DO i = 1,jpi
255          !
256          IF (Grid%bathy_meter(i,j) .EQ. 0.0 ) THEN
257             Grid%bathy_level(i,j)=0
258          ELSE
259             !
260             k1=4
261             DO WHILE (k1 .LT. (N-1))
262                IF ((Grid%bathy_meter(i,j).GE.gdepw(k1)) &
263                     .AND.(Grid%bathy_meter(i,j).LE.gdepw(k1+1))) EXIT
264                k1=k1+1
265             END DO
266             Grid%bathy_level(i,j)=k1
267             !
268          ENDIF
269          !
270       END DO
271    END DO
272    !
273  END SUBROUTINE meter_to_levels
274  !
275  !!
276  !****************************************************************
277  !   subroutine levels_to_meter             *
278  !                        *
279  ! subroutine to convert bathymetry in meters to bathymetry   *
280  ! in vertical levels                 *
281  !                        *
282  ! - input/output :                *
283  !     Grid : grid where conversion is required         *
284  !                        *
285  !various input parameters come from namelist.input files  *
286  !****************************************************************
287  !
288  SUBROUTINE levels_to_meter(Grid)
289    !
290    IMPLICIT NONE
291    !
292    REAL*8 :: za2,za1,za0,zsur,zacr,zkth,zacr2,zkth2,zmin,zmax
293    TYPE(Coordinates) :: Grid
294    INTEGER :: i,j
295    INTEGER, DIMENSION(1) :: k
296    INTEGER :: k1,ji,jj,jpi,jpj
297    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
298    !       
299    WRITE(*,*) 'convert bathymetry in meters for smoothing'
300    !
301    jpi = SIZE(Grid%bathy_level,1)
302    jpj = SIZE(Grid%bathy_level,2)
303    !             
304    !             
305    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
306         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
307       !   
308       za1=( ppdzmin - pphmax / (N-1) )          &
309            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
310            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
311            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
312
313       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
314       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
315       !
316    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
317         pa0.NE.0 .AND. pa1.NE.0 ) THEN
318       !       
319       zsur = psur
320       za0  = pa0
321       za1  = pa1
322       za2  = pa2
323       !
324    ELSE
325       !       
326       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
327       WRITE(*,*) ' '
328       WRITE(*,*) 'please check values of variables'
329       WRITE(*,*) 'in namelist vertical_grid section'
330       WRITE(*,*) ' ' 
331       STOP     
332       !       
333    ENDIF
334
335    zacr = ppacr
336    zkth = ppkth
337    zacr2 = ppacr2
338    zkth2 = ppkth2   
339    !
340    ALLOCATE(gdepw(N),gdept(N),e3w(N),e3t(N))
341    !
342    IF( ppkth == 0. ) THEN            !  uniform vertical grid
343       za1 = pphmax / FLOAT(N-1) 
344       DO i = 1, N
345          gdepw(i) = ( i - 1   ) * za1
346          gdept(i) = ( i - 0.5 ) * za1
347          e3w  (i) =  za1
348          e3t  (i) =  za1
349       END DO
350    ELSE                            ! Madec & Imbard 1996 function
351       IF( .NOT. ldbletanh ) THEN
352          DO i = 1,N
353             !
354             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
355             gdept(i) = (zsur+za0*(i+0.5)+za1*zacr*LOG(COSH(((i+0.5)-zkth)/zacr)))
356             e3w(i)   = (za0 + za1 * TANH((i-zkth)/zacr))
357             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
358             !
359          END DO
360       ELSE
361          DO i = 1,N
362             ! Double tanh function
363             gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
364                &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
365             gdept(i) = ( zsur + za0*(i+0.5) + za1 * zacr * LOG ( COSH( ((i+0.5)-zkth ) / zacr  ) )    &
366                &                            + za2 * zacr2* LOG ( COSH( ((i+0.5)-zkth2) / zacr2 ) )  )
367             e3w  (i) =          za0         + za1        * TANH(       (i-zkth ) / zacr  )            &
368                &                            + za2        * TANH(       (i-zkth2) / zacr2 )
369             e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
370                &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
371         END DO
372       ENDIF
373    ENDIF
374    !
375    gdepw(1) = 0.0 
376    !
377    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
378       !
379       DO i = 1, N-1
380          e3t(i) = gdepw(i+1)-gdepw(i)
381       END DO
382       e3t(N) = e3t(N-1)
383
384       DO i = 2, N
385          e3w(i) = gdept(i) - gdept(i-1)
386       END DO
387       e3w(1  ) = 2. * (gdept(1) - gdepw(1))
388    END IF
389    !
390    IF(.NOT. ASSOCIATED(Grid%bathy_meter)) THEN
391       ALLOCATE(Grid%bathy_meter(jpi,jpj)) 
392    ELSE
393       IF( ANY(SHAPE(Grid%bathy_meter)/=(/jpi,jpj/)) ) THEN   
394          DEALLOCATE(Grid%bathy_meter)   
395          ALLOCATE(Grid%bathy_meter(jpi,jpj))     
396       ENDIF
397    ENDIF
398    !     
399    DO jj = 1, jpj
400       DO ji= 1, jpi
401          !
402          Grid%bathy_meter(ji,jj) = gdepw( INT( Grid%bathy_level(ji,jj) ) + 1 )
403          !
404       END DO
405    END DO
406    !
407  END SUBROUTINE levels_to_meter
408  !
409
410  !****************************************************************
411  !   subroutine smooth_topo              *
412  !                        *
413  ! subroutine to smooth a given bathymetry (in meters)     *
414  !  hanning filter is used (smoothing criterion : rfactor) *
415  !                        *
416  ! - input/output :                *
417  !     h : bathymetry                 *
418  !                        *
419  !various input parameters are stored in namelist.input files *
420  !****************************************************************
421  !
422  SUBROUTINE smooth_topo(h,nbiter)
423    !
424    IMPLICIT NONE
425    !
426    REAL*8, DIMENSION(:,:) :: h
427    REAL*8 :: hmin,cff,nu,r
428    REAL*8, DIMENSION(:,:), ALLOCATABLE :: rx,ry,cx,cy,f_x,f_y
429    INTEGER :: Mm,Mmm,Lm,Lmm,M,L,nbiter,i,j
430    REAL*8,DIMENSION(:,:),ALLOCATABLE :: maskedtopo
431    !
432    M = SIZE(h,1)
433    L = SIZE(h,2)     
434    !       
435    ALLOCATE(cx(M,L),cy(M,L))
436    ALLOCATE(rx(M,L),ry(M,L))
437    ALLOCATE(f_x(M,L),f_y(M,L))
438    ALLOCATE(maskedtopo(M,L))
439    !
440    WRITE(*,*) ''
441    WRITE(*,*) 'smooth the topography (Hanning filter)'
442    WRITE(*,*) 'slope parameter = ',smoothing_factor
443    !
444    hmin = 1.1
445    WHERE(h <= hmin)
446       h = hmin
447    END WHERE
448    !       
449    WHERE (h == hmin)
450       maskedtopo = 0.
451    ELSEWHERE
452       maskedtopo = 1.
453    END WHERE
454    !
455    Mm = M-1
456    Mmm = Mm - 1
457    Lm = L-1
458    Lmm = Lm - 1
459    cff = 0.8
460    nu = 3.0/16.0
461    rx=0.
462    ry=0.     
463    CALL rfact(h,rx,ry,maskedtopo)   
464    r = MAX(MAXVAL(rx),MAXVAL(ry))
465    h = LOG(h)         
466    nbiter = 0
467    !       
468    DO WHILE (r.GT.smoothing_factor .AND. nbiter < 500  )         
469       !       
470       nbiter=nbiter+1       
471       WHERE(rx > cff*smoothing_factor)
472          cx = 1
473       ELSEWHERE
474          cx = 0
475       END WHERE
476       CALL hanningx(cx,maskedtopo)     
477       WHERE(ry > cff*smoothing_factor)
478          cy = 1
479       ELSEWHERE
480          cy = 0
481       END WHERE
482       CALL hanningy(cy,maskedtopo)     
483       CALL FX(h,f_x,cx,maskedtopo)
484       CALL FY(h,f_y,cy,maskedtopo)       
485       h(2:Mm,2:Lm) = h(2:Mm,2:Lm) + maskedtopo(2:Mm,2:Lm)*nu *                 &
486            ((f_x(2:Mm,3:L)-f_x(2:Mm,2:Lm)) + &     
487            (f_y(3:M,2:Lm)-f_y(2:Mm,2:Lm)))       
488       CALL rfact(EXP(h),rx,ry,maskedtopo)
489       r = MAX(MAXVAL(rx(2:Mm,2:L)),MAXVAL(ry(2:M,2:Lm)))
490       !
491    END DO
492    !       
493    WRITE(*,*) 'iterations = ',nbiter
494    WRITE(*,*) ''
495    h = EXP(h)               
496    WHERE( ABS(h-hmin) <= 0.001 )
497       h = 0.
498    END WHERE
499    DEALLOCATE(rx,ry,cx,cy,f_x,f_y,maskedtopo)
500    !
501  END SUBROUTINE smooth_topo
502  !       
503  !************************************************************************
504  ! subroutine hanning(bathy_meter)
505  !************************************************************************
506  !
507  SUBROUTINE hanning(h,maskedtopo)
508    !
509    IMPLICIT NONE
510    !
511    REAL*8, DIMENSION(:,:) :: h,maskedtopo
512    !
513    INTEGER :: Mm,Mmm,Lm,Lmm,M,L
514    !     
515    M = SIZE(h,1)
516    L = SIZE(h,2)
517    Mm = M-1
518    Mmm = Mm - 1
519    Lm = L-1
520    Lmm = Lm - 1
521    !     
522    h(2:Mm,2:Lm) = maskedtopo(2:Mm,2:Lm)*0.125*(   h(1:Mmm,2:Lm) + &
523         h(3:M,2:Lm)   + &
524         h(2:Mm,1:Lmm) + &
525         h(2:Mm,3:L)   + &
526         4*h(2:Mm,2:Lm))+(1.-maskedtopo(2:Mm,2:Lm))*h(2:Mm,2:Lm)
527    !
528  END SUBROUTINE hanning
529  !     
530  !************************************************************************
531  ! subroutine hanningx(bathy_meter)
532  !************************************************************************
533  !
534  SUBROUTINE hanningx(h,maskedtopo)
535    !
536    IMPLICIT NONE
537    !
538    REAL*8, DIMENSION(:,:) :: h,maskedtopo
539    REAL*8, DIMENSION(:,:), ALLOCATABLE ::  htemp 
540    !
541    INTEGER :: Mm,Mmm,Lm,Lmm,M,L
542    INTEGER :: i,j
543    !     
544    M = SIZE(h,1)
545    L = SIZE(h,2)
546    Mm = M-1
547    Mmm = Mm - 1
548    Lm = L-1
549    Lmm = Lm - 1
550
551    ALLOCATE(htemp(M,L))
552    !     
553    htemp = h     
554    DO j=3,Lm
555       DO i=2,Mm
556          IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .NE.0.) THEN
557             h(i,j)=0.125*(htemp(i-1,j)+htemp(i+1,j) &
558                  +htemp(i,j+1)+htemp(i,j-1)+4.*htemp(i,j))
559          ENDIF
560       ENDDO
561    ENDDO
562    j=2
563    DO i=2,Mm
564       IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .NE.0.) THEN
565          h(i,j)=0.25*(htemp(i+1,j)+htemp(i-1,j)+2.*htemp(i,j))
566       ENDIF
567    ENDDO
568    j=L
569    DO i=2,Mm
570       IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .NE.0.) THEN
571          h(i,j)=0.25*(htemp(i+1,j)+htemp(i-1,j)+2.*htemp(i,j))
572       ENDIF
573    ENDDO
574    DEALLOCATE(htemp)
575    !
576  END SUBROUTINE hanningx
577
578  !************************************************************************
579  ! subroutine hanning(bathy_meter)
580  !************************************************************************
581  !
582  SUBROUTINE hanningy(h,maskedtopo)
583    !
584    IMPLICIT NONE
585    !
586    REAL*8, DIMENSION(:,:) :: h,maskedtopo
587    REAL*8, DIMENSION(:,:), ALLOCATABLE ::  htemp 
588    !
589    INTEGER :: Mm,Mmm,Lm,Lmm,M,L
590    INTEGER :: i,j
591    !     
592    M = SIZE(h,1)
593    L = SIZE(h,2)
594    Mm = M-1
595    Mmm = Mm - 1
596    Lm = L-1
597    Lmm = Lm - 1     
598    ALLOCATE(htemp(M,L))
599    !     
600    htemp = h
601
602    DO j=2,Lm
603       DO i=3,Mm
604          IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .NE.0.) THEN
605             h(i,j)=0.125*(htemp(i-1,j)+htemp(i+1,j) &
606                  +htemp(i,j+1)+htemp(i,j-1)+4.*htemp(i,j))
607          ENDIF
608       ENDDO
609    ENDDO
610
611    i=2
612    DO j=2,Lm
613       IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .NE.0.) THEN
614          h(i,j)=0.25*(htemp(i,j+1)+htemp(i,j-1)+2.*htemp(i,j))
615       ENDIF
616    ENDDO
617
618    i=M
619    DO j=2,Lm
620       IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .NE.0.) THEN
621          h(i,j)=0.25*(htemp(i,j+1)+htemp(i,j-1)+2.*htemp(i,j))
622       ENDIF
623    ENDDO
624
625    DEALLOCATE(htemp)
626    !
627  END SUBROUTINE hanningy
628
629  !     
630  !************************************************************************
631  ! subroutine FX(bathy_meter,fx)
632  !************************************************************************
633  !
634  SUBROUTINE FX(h,f,c,maskedtopo)
635    !
636    IMPLICIT NONE
637    !
638    REAL*8, DIMENSION(:,:)  :: h,c
639    REAL*8, DIMENSION(:,:)  :: f,maskedtopo
640    REAL*8, DIMENSION(SIZE(h,1),SIZE(h,2))  :: floc
641    !
642    INTEGER :: Mm,Mmm,Lm,Lmm,M,L,i,j
643    !
644    f = 0.0     
645    M = SIZE(h,1)
646    L = SIZE(h,2)
647    Mm = M-1
648    Mmm = Mm - 1
649    Lm = L-1
650    Lmm = Lm - 1
651    floc = 0. 
652
653    DO j=2,L
654       DO i=1,M
655
656          IF ((maskedtopo(i,j)*maskedtopo(i,j-1)).EQ.0.) THEN
657             floc(i,j)=0.
658          ELSEIF ((i.EQ.1).OR.(i.EQ.M)) THEN     
659             floc(i,j)=(7./12.)*(h(i,j)-h(i,j-1))             
660          ELSEIF ((maskedtopo(i-1,j)*maskedtopo(i-1,j-1)).EQ.0.) THEN
661             floc(i,j)=(7./12.)*(h(i,j)-h(i,j-1))   
662          ELSEIF ((maskedtopo(i+1,j)*maskedtopo(i+1,j-1)).EQ.0.) THEN
663             floc(i,j)=(7./12.)*(h(i,j)-h(i,j-1))       
664          ELSE               
665             floc(i,j)=(5./12.)*(h(i,j)-h(i,j-1)) &
666                  +(1./12.)*(h(i-1,j)-h(i-1,j-1)+h(i+1,j)-h(i+1,j-1))
667          ENDIF
668       ENDDO
669    ENDDO
670    !     
671    DO j = 1,L
672       DO i = 1,M
673          f(i,j) = c(i,j)*floc(i,j) 
674       END DO
675    END DO
676    !                   
677    !
678  END SUBROUTINE FX
679  !
680  !
681  !     
682  !************************************************************************
683  ! subroutine FY(bathy_meter,fy)
684  !************************************************************************
685  !
686  SUBROUTINE FY(h,f,c,maskedtopo)
687    !
688    IMPLICIT NONE
689    !
690    REAL*8, DIMENSION(:,:) :: h,c
691    REAL*8, DIMENSION(:,:) :: f,maskedtopo
692    REAL*8, DIMENSION(SIZE(h,1),SIZE(h,2))  :: floc
693    INTEGER :: Mm,Mmm,Lm,Lmm,M,L,i,j
694    f=0.0 
695    !           
696    M = SIZE(h,1)
697    L = SIZE(h,2)
698    Mm = M-1
699    Mmm = Mm - 1
700    Lm = L-1
701    Lmm = Lm - 1
702    !
703    floc = 0.
704
705    DO j=1,L
706       DO i=2,M
707          IF ((maskedtopo(i,j)*maskedtopo(i-1,j)).EQ.0.) THEN
708             floc(i,j) = 0.
709          ELSEIF ((j.EQ.1).OR.(j.EQ.L)) THEN
710             floc(i,j)=(7./12.)*(h(i,j)-h(i-1,j))
711          ELSEIF ((maskedtopo(i,j-1)*maskedtopo(i-1,j-1)).EQ.0.) THEN
712             floc(i,j)=(7./12.)*(h(i,j)-h(i-1,j))
713          ELSEIF ((maskedtopo(i,j+1)*maskedtopo(i-1,j+1)).EQ.0.) THEN
714             floc(i,j)=(7./12.)*(h(i,j)-h(i-1,j))
715          ELSE     
716             floc(i,j)=(5./12.)*(h(i,j)-h(i-1,j)) &
717                  +(1./12.)*(h(i,j-1)-h(i-1,j-1)+h(i,j+1)-h(i-1,j+1))
718          ENDIF
719       ENDDO
720    ENDDO
721    !
722    DO j = 1,L
723       DO i = 1,M
724          f(i,j) = c(i,j)*floc(i,j) 
725       END DO
726    END DO
727    !     
728  END SUBROUTINE FY
729  !
730  !
731  !****************************************************************
732  !   subroutine rfact                 *
733  !                        *
734  ! subroutine to check if smoothing criterion        *
735  !                     is verified everywhere        *
736  !                        *
737  ! - input :                    *
738  !     h : bathymetry                 *
739  ! - ouput :                    *
740  !    rx,ry : delta(theta)/theta in x and y directions     *
741  !****************************************************************
742  !
743  SUBROUTINE rfact(h,rx,ry,maskedtopo)
744    !
745    IMPLICIT NONE
746    !
747    REAL*8, DIMENSION(:,:)  :: h
748    REAL*8, DIMENSION(:,:)  :: rx,ry
749    REAL*8, DIMENSION(:,:)  :: maskedtopo
750    INTEGER M,L,i,j,Mm,Mmm,Lm,Lmm
751    !
752    M = SIZE(h,1)
753    L = SIZE(h,2)
754    Mm = M-1
755    Mmm = Mm - 1
756    Lm = L-1
757    Lmm = Lm - 1
758    !
759    rx=0.
760    ry=0.
761    !     
762    DO j=2,L
763       DO i=1,M       
764          rx(i,j) = ABS(h(i,j)-h(i,j-1))/(h(i,j)+h(i,j-1))
765          IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .EQ.0.) THEN
766             rx(i,j)=0.
767          ENDIF
768       ENDDO
769    ENDDO
770    !
771    DO j=1,L
772       DO i=2,M
773          ry(i,j) = ABS(h(i,j)-h(i-1,j))/(h(i,j)+h(i-1,j))
774          IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .EQ.0.) THEN
775             ry(i,j)=0.
776          ENDIF
777       ENDDO
778    ENDDO
779    !
780  END SUBROUTINE rfact
781  !
782  !
783  !****************************************************************
784  !   subroutine Update_Parent_Bathy            *
785  !                        *
786  ! (if desired) subroutine to update parent grid bathymetry   *
787  ! for consistency with fine grid bathymetry         *
788  !                        *
789  ! if a given coarse grid point is masked and one of the     *
790  ! child grid points contained in this coarse cell is not masked *
791  ! the corresponding coarse grid point is unmasked with rn_hmin  *
792  ! value                            *
793  !                        *
794  ! - input :                    *
795  !     G0,G1 : both grids involved          *
796  ! - ouput :                    *
797  !    G0 parent grid containing updated bathymetry      *
798  !****************************************************************
799  !
800  !
801  SUBROUTINE Update_Parent_Bathy( G0,G1 )
802    !
803    IMPLICIT NONE
804
805    TYPE(coordinates) :: G0,G1
806    INTEGER :: ii,jj,jk,ipt,jpt,diff,indx,indy,bornex,borney,bornex2,borney2
807    !
808    REAL*8 :: za2,za1,za0,zsur,zacr,zkth,zacr2,zkth2,zmin
809    INTEGER :: i,j
810    INTEGER :: k1
811    INTEGER :: compt, compt_oce
812    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
813    !                     
814    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
815         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
816       !   
817       za1=( ppdzmin - pphmax / (N-1) )          &
818            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
819            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
820            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
821
822       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
823       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
824       !
825    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
826         pa0.NE.0 .AND. pa1.NE.0 ) THEN
827       !       
828       zsur = psur
829       za0  = pa0
830       za1  = pa1
831       za2  = pa2
832       !
833    ELSE
834       !       
835       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
836       WRITE(*,*) ' '
837       WRITE(*,*) 'please check values of variables'
838       WRITE(*,*) 'in namelist vertical_grid section'
839       WRITE(*,*) ' ' 
840       STOP     
841       !       
842    ENDIF
843
844    zacr = ppacr
845    zkth = ppkth
846    zacr2 = ppacr2
847    zkth2 = ppkth2   
848    !
849    ALLOCATE(gdepw(N),gdept(N),e3w(N),e3t(N))
850    !
851    IF( ppkth == 0. ) THEN            !  uniform vertical grid
852       za1 = pphmax / FLOAT(N-1) 
853       DO i = 1, N
854          gdepw(i) = ( i - 1   ) * za1
855          gdept(i) = ( i - 0.5 ) * za1
856          e3w  (i) =  za1
857          e3t  (i) =  za1
858       END DO
859    ELSE                            ! Madec & Imbard 1996 function
860       IF( .NOT. ldbletanh ) THEN
861          DO i = 1,N
862             !
863             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
864             gdept(i) = (zsur+za0*(i+0.5)+za1*zacr*LOG(COSH(((i+0.5)-zkth)/zacr)))
865             e3w(i)   = (za0 + za1 * TANH((i-zkth)/zacr))
866             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
867             !
868          END DO
869       ELSE
870          DO i = 1,N
871             ! Double tanh function
872             gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
873                &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
874             gdept(i) = ( zsur + za0*(i+0.5) + za1 * zacr * LOG ( COSH( ((i+0.5)-zkth ) / zacr  ) )    &
875                &                            + za2 * zacr2* LOG ( COSH( ((i+0.5)-zkth2) / zacr2 ) )  )
876             e3w  (i) =          za0         + za1        * TANH(       (i-zkth ) / zacr  )            &
877                &                            + za2        * TANH(       (i-zkth2) / zacr2 )
878             e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
879                &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
880          END DO
881       ENDIF
882    ENDIF
883    !
884    gdepw(1)=0. 
885    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
886       !
887       DO i = 1, N-1
888          e3t(i) = gdepw(i+1)-gdepw(i)
889       END DO
890       e3t(N) = e3t(N-1)
891
892       DO i = 2, N
893          e3w(i) = gdept(i) - gdept(i-1)
894       END DO
895       e3w(1  ) = 2. * (gdept(1) - gdepw(1))
896    END IF
897    !
898    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
899    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
900    ENDIF
901    zmin = gdepw(i+1)
902    !     
903    diff = 0
904    IF(MOD(rho,2) .EQ. 0) diff = 1
905    !
906    indx = 1+nbghostcellsfine + CEILING(irafx/2.0) + diff     
907    indy = 1+nbghostcellsfine + CEILING(irafy/2.0) + diff
908    bornex = 1+nbghostcellsfine + CEILING(irafx/2.0) + diff
909    borney = 1+nbghostcellsfine + CEILING(irafy/2.0) + diff
910    bornex2 = nxfin - (nbghostcellsfine) - CEILING(irafx/2.0) 
911    borney2 = nyfin - (nbghostcellsfine) - CEILING(irafy/2.0) 
912    !
913    DO j = borney,borney2,irafy
914
915       jpt = jmin + 1 + nbghostcellsfine + (j-indy)/irafy
916       IF(j<=1) jpt = jmin + 1
917
918       DO i = bornex,bornex2,irafx
919
920          ipt = imin + 1 + nbghostcellsfine + (i-indx)/irafx
921          IF(i<=1) ipt = imin + 1
922          IF ((i==bornex).AND.(j==borney)) print *, 'Coarse grid Corner', ipt,jpt
923          !   
924          G0%Bathy_meter(ipt,jpt) = 0.
925          compt = 0 
926          compt_oce = 0
927          DO jj = j-FLOOR(irafy/2.0),j+FLOOR(irafy/2.0)-diff
928             DO ii = i-FLOOR(irafx/2.0),i+FLOOR(irafx/2.0)-diff
929                G0%Bathy_meter(ipt,jpt) = G0%Bathy_meter(ipt,jpt) + G1%Bathy_meter(ii,jj)
930                IF (G1%Bathy_meter(ii,jj)>0.) compt_oce = compt_oce + 1
931                compt = compt + 1
932             END DO
933          END DO
934!          IF (compt_oce > 0.5*irafx*irafy)    G0%Bathy_meter(ipt,jpt) = G0%Bathy_meter(ipt,jpt) / FLOAT(compt)
935          G0%Bathy_meter(ipt,jpt) = G0%Bathy_meter(ipt,jpt) / FLOAT(compt)
936          IF(G0%Bathy_meter(ipt,jpt).GT.0.)   G0%Bathy_meter(ipt,jpt) = MAX(G0%Bathy_meter(ipt,jpt),zmin)
937          !
938          G0%wgt(ipt,jpt) = 1. ! Flag for output
939       END DO
940    END DO
941
942    !
943    WRITE(*,*) ' Number of coarse grid points updated = ',compt
944    WRITE(*,*) '---------------------------------'
945    !     
946  END SUBROUTINE Update_Parent_Bathy
947
948END MODULE agrif_connect_topo
Note: See TracBrowser for help on using the repository browser.