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/UKMO/dev_r5518_rm_um_cpl/NEMOGCM/TOOLS/NESTING/src – NEMO

source: branches/UKMO/dev_r5518_rm_um_cpl/NEMOGCM/TOOLS/NESTING/src/agrif_connect_topo.f90 @ 5855

Last change on this file since 5855 was 5855, checked in by jcastill, 8 years ago

Clearing SVN keywords

File size: 21.9 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 :: za1,za0,zsur,zacr,zkth,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       !
137    ELSE
138       !       
139       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
140       WRITE(*,*) ' '
141       WRITE(*,*) 'please check values of variables'
142       WRITE(*,*) 'in namelist vertical_grid section'
143       WRITE(*,*) ' ' 
144       STOP     
145       !       
146    ENDIF
147
148    zacr = ppacr
149    zkth = ppkth 
150
151    !
152    ALLOCATE(gdepw(N),gdept(N),e3w(N),e3t(N))
153    !
154    DO i = 1,N
155       gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
156       gdept(i) = (zsur+za0*(i+0.5)+za1*zacr*LOG(COSH(((i+0.5)-zkth)/zacr)))
157       e3w(i)   = (za0 + za1 * TANH((i-zkth)/zacr))
158       e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
159    END DO
160    !
161    gdepw(1) = 0.0
162    zmax = gdepw(N) + e3t(N)
163    zmin = gdepw(4)
164    !
165    IF ( .NOT. ASSOCIATED(Grid%bathy_level)) &
166         ALLOCATE(Grid%bathy_level(jpi,jpj))   
167    !
168    Grid%bathy_level = N-1
169    !
170    DO jj = 1, jpj
171       DO ji= 1, jpi
172          IF( Grid%bathy_meter(ji,jj) <= 0. )   &
173               Grid%bathy_level(ji,jj) = INT( Grid%bathy_meter(ji,jj) )
174       END DO
175    END DO
176    !
177    DO jj = 1, jpj
178       DO ji= 1, jpi
179          IF( Grid%bathy_meter(ji,jj) <= 0. ) THEN
180             Grid%bathy_meter(ji,jj) = 0.e0
181          ELSE
182             Grid%bathy_meter(ji,jj) = MAX( Grid%bathy_meter(ji,jj), zmin )
183             Grid%bathy_meter(ji,jj) = MIN( Grid%bathy_meter(ji,jj), zmax )
184          ENDIF
185       END DO
186    END DO
187    !
188    !
189    !
190    DO j = 1,jpj
191       DO i = 1,jpi
192          !
193          IF (Grid%bathy_meter(i,j) .EQ. 0.0 ) THEN
194             Grid%bathy_level(i,j)=0.
195          ELSE
196             !
197             k1=4
198             DO WHILE (k1 .LT. (N-1))
199                IF ((Grid%bathy_meter(i,j).GE.gdepw(k1)) &
200                     .AND.(Grid%bathy_meter(i,j).LE.gdepw(k1+1))) EXIT
201                k1=k1+1
202             END DO
203             Grid%bathy_level(i,j)=k1
204             !
205          ENDIF
206          !
207       END DO
208    END DO
209    !
210  END SUBROUTINE meter_to_levels
211  !
212  !!
213  !****************************************************************
214  !   subroutine levels_to_meter             *
215  !                        *
216  ! subroutine to convert bathymetry in meters to bathymetry   *
217  ! in vertical levels                 *
218  !                        *
219  ! - input/output :                *
220  !     Grid : grid where conversion is required         *
221  !                        *
222  !various input parameters come from namelist.input files  *
223  !****************************************************************
224  !
225  SUBROUTINE levels_to_meter(Grid)
226    !
227    IMPLICIT NONE
228    !
229    REAL*8 :: za1,za0,zsur,zacr,zkth,zmin,zmax
230    TYPE(Coordinates) :: Grid
231    INTEGER :: i,j
232    INTEGER, DIMENSION(1) :: k
233    INTEGER :: k1,ji,jj,jpi,jpj
234    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
235    !       
236    WRITE(*,*) 'convert bathymetry in meters for smoothing'
237    !
238    jpi = SIZE(Grid%bathy_level,1)
239    jpj = SIZE(Grid%bathy_level,2)
240    !             
241    !             
242    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
243         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
244       !   
245       za1=( ppdzmin - pphmax / (N-1) )          &
246            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
247            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
248            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
249
250       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
251       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
252       !
253    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
254         pa0.NE.0 .AND. pa1.NE.0 ) THEN
255       !       
256       zsur = psur
257       za0  = pa0
258       za1  = pa1
259       !
260    ELSE
261       !       
262       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
263       WRITE(*,*) ' '
264       WRITE(*,*) 'please check values of variables'
265       WRITE(*,*) 'in namelist vertical_grid section'
266       WRITE(*,*) ' '
267       STOP     
268       !       
269    ENDIF
270    !
271    zacr = ppacr
272    zkth = ppkth 
273
274    !
275    ALLOCATE(gdepw(N),gdept(N),e3w(N),e3t(N))
276    !
277    DO i = 1,N
278       !
279       gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
280       gdept(i) = (zsur+za0*(i+0.5)+za1*zacr*LOG(COSH(((i+0.5)-zkth)/zacr)))
281       e3w(i)   = (za0 + za1 * TANH((i-zkth)/zacr))
282       e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
283    END DO
284    !
285    gdepw(1) = 0.0 
286    !
287    IF(.NOT. ASSOCIATED(Grid%bathy_meter)) THEN
288       ALLOCATE(Grid%bathy_meter(jpi,jpj)) 
289    ELSE
290       IF( ANY(SHAPE(Grid%bathy_meter)/=(/jpi,jpj/)) ) THEN   
291          DEALLOCATE(Grid%bathy_meter)   
292          ALLOCATE(Grid%bathy_meter(jpi,jpj))     
293       ENDIF
294    ENDIF
295    !     
296    DO jj = 1, jpj
297       DO ji= 1, jpi
298          !
299          Grid%bathy_meter(ji,jj) = gdepw( INT( Grid%bathy_level(ji,jj) ) + 1 )
300          !
301       END DO
302    END DO
303    !
304  END SUBROUTINE levels_to_meter
305  !
306
307  !****************************************************************
308  !   subroutine smooth_topo              *
309  !                        *
310  ! subroutine to smooth a given bathymetry (in meters)     *
311  !  hanning filter is used (smoothing criterion : rfactor) *
312  !                        *
313  ! - input/output :                *
314  !     h : bathymetry                 *
315  !                        *
316  !various input parameters are stored in namelist.input files *
317  !****************************************************************
318  !
319  SUBROUTINE smooth_topo(h,nbiter)
320    !
321    IMPLICIT NONE
322    !
323    REAL*8, DIMENSION(:,:) :: h
324    REAL*8 :: hmin,cff,nu,r
325    REAL*8, DIMENSION(:,:), ALLOCATABLE :: rx,ry,cx,cy,f_x,f_y
326    INTEGER :: Mm,Mmm,Lm,Lmm,M,L,nbiter,i,j
327    REAL*8,DIMENSION(:,:),ALLOCATABLE :: maskedtopo
328    !
329    M = SIZE(h,1)
330    L = SIZE(h,2)     
331    !       
332    ALLOCATE(cx(M,L),cy(M,L))
333    ALLOCATE(rx(M,L),ry(M,L))
334    ALLOCATE(f_x(M,L),f_y(M,L))
335    ALLOCATE(maskedtopo(M,L))
336    !
337    WRITE(*,*) ''
338    WRITE(*,*) 'smooth the topography (Hanning filter)'
339    WRITE(*,*) 'slope parameter = ',smoothing_factor
340    !
341    hmin = 1.1
342    WHERE(h <= hmin)
343       h = hmin
344    END WHERE
345    !       
346    WHERE (h == hmin)
347       maskedtopo = 0.
348    ELSEWHERE
349       maskedtopo = 1.
350    END WHERE
351    !
352    Mm = M-1
353    Mmm = Mm - 1
354    Lm = L-1
355    Lmm = Lm - 1
356    cff = 0.8
357    nu = 3.0/16.0
358    rx=0.
359    ry=0.     
360    CALL rfact(h,rx,ry,maskedtopo)   
361    r = MAX(MAXVAL(rx),MAXVAL(ry))
362    h = LOG(h)         
363    nbiter = 0
364    !       
365    DO WHILE (r.GT.smoothing_factor .AND. nbiter < 500  )         
366       !       
367       nbiter=nbiter+1       
368       WHERE(rx > cff*smoothing_factor)
369          cx = 1
370       ELSEWHERE
371          cx = 0
372       END WHERE
373       CALL hanningx(cx,maskedtopo)     
374       WHERE(ry > cff*smoothing_factor)
375          cy = 1
376       ELSEWHERE
377          cy = 0
378       END WHERE
379       CALL hanningy(cy,maskedtopo)     
380       CALL FX(h,f_x,cx,maskedtopo)
381       CALL FY(h,f_y,cy,maskedtopo)       
382       h(2:Mm,2:Lm) = h(2:Mm,2:Lm) + maskedtopo(2:Mm,2:Lm)*nu *                 &
383            ((f_x(2:Mm,3:L)-f_x(2:Mm,2:Lm)) + &     
384            (f_y(3:M,2:Lm)-f_y(2:Mm,2:Lm)))       
385       CALL rfact(EXP(h),rx,ry,maskedtopo)
386       r = MAX(MAXVAL(rx(2:Mm,2:L)),MAXVAL(ry(2:M,2:Lm)))
387       !
388    END DO
389    !       
390    WRITE(*,*) 'iterations = ',nbiter
391    WRITE(*,*) ''
392    h = EXP(h)               
393    WHERE( ABS(h-hmin) <= 0.001 )
394       h = 0.
395    END WHERE
396    DEALLOCATE(rx,ry,cx,cy,f_x,f_y,maskedtopo)
397    !
398  END SUBROUTINE smooth_topo
399  !       
400  !************************************************************************
401  ! subroutine hanning(bathy_meter)
402  !************************************************************************
403  !
404  SUBROUTINE hanning(h,maskedtopo)
405    !
406    IMPLICIT NONE
407    !
408    REAL*8, DIMENSION(:,:) :: h,maskedtopo
409    !
410    INTEGER :: Mm,Mmm,Lm,Lmm,M,L
411    !     
412    M = SIZE(h,1)
413    L = SIZE(h,2)
414    Mm = M-1
415    Mmm = Mm - 1
416    Lm = L-1
417    Lmm = Lm - 1
418    !     
419    h(2:Mm,2:Lm) = maskedtopo(2:Mm,2:Lm)*0.125*(   h(1:Mmm,2:Lm) + &
420         h(3:M,2:Lm)   + &
421         h(2:Mm,1:Lmm) + &
422         h(2:Mm,3:L)   + &
423         4*h(2:Mm,2:Lm))+(1.-maskedtopo(2:Mm,2:Lm))*h(2:Mm,2:Lm)
424    !
425  END SUBROUTINE hanning
426  !     
427  !************************************************************************
428  ! subroutine hanningx(bathy_meter)
429  !************************************************************************
430  !
431  SUBROUTINE hanningx(h,maskedtopo)
432    !
433    IMPLICIT NONE
434    !
435    REAL*8, DIMENSION(:,:) :: h,maskedtopo
436    REAL*8, DIMENSION(:,:), ALLOCATABLE ::  htemp 
437    !
438    INTEGER :: Mm,Mmm,Lm,Lmm,M,L
439    INTEGER :: i,j
440    !     
441    M = SIZE(h,1)
442    L = SIZE(h,2)
443    Mm = M-1
444    Mmm = Mm - 1
445    Lm = L-1
446    Lmm = Lm - 1
447
448    ALLOCATE(htemp(M,L))
449    !     
450    htemp = h     
451    DO j=3,Lm
452       DO i=2,Mm
453          IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .NE.0.) THEN
454             h(i,j)=0.125*(htemp(i-1,j)+htemp(i+1,j) &
455                  +htemp(i,j+1)+htemp(i,j-1)+4.*htemp(i,j))
456          ENDIF
457       ENDDO
458    ENDDO
459    j=2
460    DO i=2,Mm
461       IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .NE.0.) THEN
462          h(i,j)=0.25*(htemp(i+1,j)+htemp(i-1,j)+2.*htemp(i,j))
463       ENDIF
464    ENDDO
465    j=L
466    DO i=2,Mm
467       IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .NE.0.) THEN
468          h(i,j)=0.25*(htemp(i+1,j)+htemp(i-1,j)+2.*htemp(i,j))
469       ENDIF
470    ENDDO
471    DEALLOCATE(htemp)
472    !
473  END SUBROUTINE hanningx
474
475  !************************************************************************
476  ! subroutine hanning(bathy_meter)
477  !************************************************************************
478  !
479  SUBROUTINE hanningy(h,maskedtopo)
480    !
481    IMPLICIT NONE
482    !
483    REAL*8, DIMENSION(:,:) :: h,maskedtopo
484    REAL*8, DIMENSION(:,:), ALLOCATABLE ::  htemp 
485    !
486    INTEGER :: Mm,Mmm,Lm,Lmm,M,L
487    INTEGER :: i,j
488    !     
489    M = SIZE(h,1)
490    L = SIZE(h,2)
491    Mm = M-1
492    Mmm = Mm - 1
493    Lm = L-1
494    Lmm = Lm - 1     
495    ALLOCATE(htemp(M,L))
496    !     
497    htemp = h
498
499    DO j=2,Lm
500       DO i=3,Mm
501          IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .NE.0.) THEN
502             h(i,j)=0.125*(htemp(i-1,j)+htemp(i+1,j) &
503                  +htemp(i,j+1)+htemp(i,j-1)+4.*htemp(i,j))
504          ENDIF
505       ENDDO
506    ENDDO
507
508    i=2
509    DO j=2,Lm
510       IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .NE.0.) THEN
511          h(i,j)=0.25*(htemp(i,j+1)+htemp(i,j-1)+2.*htemp(i,j))
512       ENDIF
513    ENDDO
514
515    i=M
516    DO j=2,Lm
517       IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .NE.0.) THEN
518          h(i,j)=0.25*(htemp(i,j+1)+htemp(i,j-1)+2.*htemp(i,j))
519       ENDIF
520    ENDDO
521
522    DEALLOCATE(htemp)
523    !
524  END SUBROUTINE hanningy
525
526  !     
527  !************************************************************************
528  ! subroutine FX(bathy_meter,fx)
529  !************************************************************************
530  !
531  SUBROUTINE FX(h,f,c,maskedtopo)
532    !
533    IMPLICIT NONE
534    !
535    REAL*8, DIMENSION(:,:)  :: h,c
536    REAL*8, DIMENSION(:,:)  :: f,maskedtopo
537    REAL*8, DIMENSION(SIZE(h,1),SIZE(h,2))  :: floc
538    !
539    INTEGER :: Mm,Mmm,Lm,Lmm,M,L,i,j
540    !
541    f = 0.0     
542    M = SIZE(h,1)
543    L = SIZE(h,2)
544    Mm = M-1
545    Mmm = Mm - 1
546    Lm = L-1
547    Lmm = Lm - 1
548    floc = 0. 
549
550    DO j=2,L
551       DO i=1,M
552
553          IF ((maskedtopo(i,j)*maskedtopo(i,j-1)).EQ.0.) THEN
554             floc(i,j)=0.
555          ELSEIF ((i.EQ.1).OR.(i.EQ.M)) THEN     
556             floc(i,j)=(7./12.)*(h(i,j)-h(i,j-1))             
557          ELSEIF ((maskedtopo(i-1,j)*maskedtopo(i-1,j-1)).EQ.0.) THEN
558             floc(i,j)=(7./12.)*(h(i,j)-h(i,j-1))   
559          ELSEIF ((maskedtopo(i+1,j)*maskedtopo(i+1,j-1)).EQ.0.) THEN
560             floc(i,j)=(7./12.)*(h(i,j)-h(i,j-1))       
561          ELSE               
562             floc(i,j)=(5./12.)*(h(i,j)-h(i,j-1)) &
563                  +(1./12.)*(h(i-1,j)-h(i-1,j-1)+h(i+1,j)-h(i+1,j-1))
564          ENDIF
565       ENDDO
566    ENDDO
567    !     
568    DO j = 1,L
569       DO i = 1,M
570          f(i,j) = c(i,j)*floc(i,j) 
571       END DO
572    END DO
573    !                   
574    !
575  END SUBROUTINE FX
576  !
577  !
578  !     
579  !************************************************************************
580  ! subroutine FY(bathy_meter,fy)
581  !************************************************************************
582  !
583  SUBROUTINE FY(h,f,c,maskedtopo)
584    !
585    IMPLICIT NONE
586    !
587    REAL*8, DIMENSION(:,:) :: h,c
588    REAL*8, DIMENSION(:,:) :: f,maskedtopo
589    REAL*8, DIMENSION(SIZE(h,1),SIZE(h,2))  :: floc
590    INTEGER :: Mm,Mmm,Lm,Lmm,M,L,i,j
591    f=0.0 
592    !           
593    M = SIZE(h,1)
594    L = SIZE(h,2)
595    Mm = M-1
596    Mmm = Mm - 1
597    Lm = L-1
598    Lmm = Lm - 1
599    !
600    floc = 0.
601
602    DO j=1,L
603       DO i=2,M
604          IF ((maskedtopo(i,j)*maskedtopo(i-1,j)).EQ.0.) THEN
605             floc(i,j) = 0.
606          ELSEIF ((j.EQ.1).OR.(j.EQ.L)) THEN
607             floc(i,j)=(7./12.)*(h(i,j)-h(i-1,j))
608          ELSEIF ((maskedtopo(i,j-1)*maskedtopo(i-1,j-1)).EQ.0.) THEN
609             floc(i,j)=(7./12.)*(h(i,j)-h(i-1,j))
610          ELSEIF ((maskedtopo(i,j+1)*maskedtopo(i-1,j+1)).EQ.0.) THEN
611             floc(i,j)=(7./12.)*(h(i,j)-h(i-1,j))
612          ELSE     
613             floc(i,j)=(5./12.)*(h(i,j)-h(i-1,j)) &
614                  +(1./12.)*(h(i,j-1)-h(i-1,j-1)+h(i,j+1)-h(i-1,j+1))
615          ENDIF
616       ENDDO
617    ENDDO
618    !
619    DO j = 1,L
620       DO i = 1,M
621          f(i,j) = c(i,j)*floc(i,j) 
622       END DO
623    END DO
624    !     
625  END SUBROUTINE FY
626  !
627  !
628  !****************************************************************
629  !   subroutine rfact                 *
630  !                        *
631  ! subroutine to check if smoothing criterion        *
632  !                     is verified everywhere        *
633  !                        *
634  ! - input :                    *
635  !     h : bathymetry                 *
636  ! - ouput :                    *
637  !    rx,ry : delta(theta)/theta in x and y directions     *
638  !****************************************************************
639  !
640  SUBROUTINE rfact(h,rx,ry,maskedtopo)
641    !
642    IMPLICIT NONE
643    !
644    REAL*8, DIMENSION(:,:)  :: h
645    REAL*8, DIMENSION(:,:)  :: rx,ry
646    REAL*8, DIMENSION(:,:)  :: maskedtopo
647    INTEGER M,L,i,j,Mm,Mmm,Lm,Lmm
648    !
649    M = SIZE(h,1)
650    L = SIZE(h,2)
651    Mm = M-1
652    Mmm = Mm - 1
653    Lm = L-1
654    Lmm = Lm - 1
655    !
656    rx=0.
657    ry=0.
658    !     
659    DO j=2,L
660       DO i=1,M       
661          rx(i,j) = ABS(h(i,j)-h(i,j-1))/(h(i,j)+h(i,j-1))
662          IF ((maskedtopo(i,j)*maskedtopo(i,j-1)) .EQ.0.) THEN
663             rx(i,j)=0.
664          ENDIF
665       ENDDO
666    ENDDO
667    !
668    DO j=1,L
669       DO i=2,M
670          ry(i,j) = ABS(h(i,j)-h(i-1,j))/(h(i,j)+h(i-1,j))
671          IF ((maskedtopo(i,j)*maskedtopo(i-1,j)) .EQ.0.) THEN
672             ry(i,j)=0.
673          ENDIF
674       ENDDO
675    ENDDO
676    !
677  END SUBROUTINE rfact
678  !
679  !
680  !****************************************************************
681  !   subroutine Update_Parent_Bathy            *
682  !                        *
683  ! (if desired) subroutine to update parent grid bathymetry   *
684  ! for consistency with fine grid bathymetry         *
685  !                        *
686  ! if a given coarse grid point is masked and one of the      *
687  ! child grid points contained in this coarse cell is not masked *
688  ! the corresponding coarse grid point is unmasked with gdepw(4) *
689  ! value                        *
690  !                        *
691  ! - input :                    *
692  !     G0,G1 : both grids involved          *
693  ! - ouput :                    *
694  !    G0 parent grid containing updated bathymetry      *
695  !****************************************************************
696  !
697  !
698  SUBROUTINE Update_Parent_Bathy( G0,G1 )
699    !
700    IMPLICIT NONE
701
702    TYPE(coordinates) :: G0,G1
703    INTEGER :: ji,jj,jk,ipt,jpt,diff,indx,indy,bornex,borney,bornex2,borney2
704    !
705    INTEGER :: ideb,jdeb,ifin,jfin
706    REAL*8 :: za1,za0,zsur,zacr,zkth,zmin
707    INTEGER :: i,j
708    INTEGER :: k1
709    INTEGER :: compt
710    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
711    !                     
712    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
713         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
714       !   
715       za1=( ppdzmin - pphmax / (N-1) )          &
716            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
717            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
718            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
719
720       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
721       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
722       !
723    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
724         pa0.NE.0 .AND. pa1.NE.0 ) THEN
725       !       
726       zsur = psur
727       za0  = pa0
728       za1  = pa1
729       !
730    ELSE
731       !       
732       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
733       WRITE(*,*) ' '
734       WRITE(*,*) 'please check values of variables'
735       WRITE(*,*) 'in namelist vertical_grid section'
736       WRITE(*,*) ' ' 
737       STOP     
738       !       
739    ENDIF
740
741    zacr = ppacr
742    zkth = ppkth 
743
744    !
745    ALLOCATE(gdepw(N),gdept(N),e3w(N),e3t(N))
746    !
747    DO i = 1,N
748       !
749       gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
750       gdept(i) = (zsur+za0*(i+0.5)+za1*zacr*LOG(COSH(((i+0.5)-zkth)/zacr)))
751       e3w(i)   = (za0 + za1 * TANH((i-zkth)/zacr))
752       e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
753    END DO
754    !
755    zmin = gdepw(4)
756    !     
757    diff = 0
758    IF(MOD(rho,2) .EQ. 0) diff = 1
759
760    ideb = 3 + (irafx - 1)/2
761    jdeb = 3 + (irafy - 1)/2
762
763    jfin = nyfin - 2 - (irafy)/2
764    ifin = nxfin - 2 - (irafx)/2
765    !
766    WRITE(*,*) '---------------------------------'
767    WRITE(*,*) 'Parent grid bathymetry update ...'
768    compt = 0
769    !
770    DO jj = jdeb,jfin,irafy
771       jpt = jmin + 3 - 1 + (jj - jdeb) / irafy
772       DO ji = ideb,ifin,irafx
773          ipt = imin + 3 - 1 + (ji - ideb) / irafx
774          !
775          IF( G0%Bathy_meter(ipt,jpt) .LE. 0. .AND.                      &
776               MAXVAL(G1%Bathy_meter(ji-irafx/2+diff:ji+irafx/2, &
777               jj-irafy/2+diff:jj+irafy/2)).GT. 0.  ) THEN
778             G0%Bathy_meter(ipt,jpt) = zmin
779             !
780             compt = compt + 1
781             !
782          ENDIF
783          !
784       ENDDO
785    ENDDO
786    !
787    WRITE(*,*) ' Number of coarse grid points updated = ',compt
788    WRITE(*,*) '---------------------------------'
789    !     
790  END SUBROUTINE Update_Parent_Bathy
791
792END MODULE agrif_connect_topo
Note: See TracBrowser for help on using the repository browser.