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 branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/TOOLS/NESTING/src/agrif_connect_topo.f90 @ 9149

Last change on this file since 9149 was 9149, checked in by jchanut, 7 years ago

NESTING TOOLS
Add new e3t definition, #1981
Correct transition weights, #1998
Correct arithmetic and median average interpolations, #1999
Changed number of coarse grid bathymetry points inside child grid domain (2 now, instead of 3)

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