New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_connect_topo.f90 in utils/tools/NESTING/src – NEMO

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

Last change on this file since 12264 was 12264, checked in by clem, 5 years ago

finishing cleaning the nesting tools (as far as I could). I will stop there
  • Property svn:keywords set to Id
File size: 28.8 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 jj = 1,jpj
254       DO ji = 1,jpi
255          !
256          IF (Grid%bathy_meter(ji,jj) .EQ. 0.0 ) THEN
257             Grid%bathy_level(ji,jj)=0
258          ELSE
259             !
260             k1=2  ! clem: minimum levels = 4 ???
261             DO WHILE (k1 .LT. (N-1))
262                IF ((Grid%bathy_meter(ji,jj).GE.gdepw(k1)) &
263                     .AND.(Grid%bathy_meter(ji,jj).LE.gdepw(k1+1))) EXIT
264                k1=k1+1
265             END DO
266             Grid%bathy_level(ji,jj)=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.