source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/TOOLS/NESTING/src/agrif_connect_topo.f90 @ 6204

Last change on this file since 6204 was 6204, checked in by cetlod, 6 years ago

back the nemo_v3_6_STABLE_XIOS2 branch into 3_6_STABLE, including bugfixes, XIOS2 and new AGRIF

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