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_partial_steps.f90 in NEMO/trunk/NEMOGCM/TOOLS/NESTING/src – NEMO

source: NEMO/trunk/NEMOGCM/TOOLS/NESTING/src/agrif_partial_steps.f90 @ 9594

Last change on this file since 9594 was 9166, checked in by jchanut, 6 years ago

NESTING TOOLS:
Fixes to account for user defined number of ghostcells - still set to 1 to be consistent with NEMO
Set child grid bathymetry near boundaries to nearest neighbor interpolation from parent
Update Coarse grid bathymetry so that each cell volume matches child grid average

  • Property svn:keywords set to Id
File size: 28.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!
8MODULE agrif_partial_steps
9  !
10  USE agrif_types
11CONTAINS
12
13
14
15
16
17  !
18  !************************************************************************
19  !                                                            *
20  ! MODULE  AGRIF_PARTIAL_STEPS                                      *
21  !                           *
22  !************************************************************************
23
24
25  !************************************************************************
26  !                           *
27  ! Subroutine get_partial_steps                *
28  !                           *
29  ! subroutine to compute gdepw_ps on the input grid (based on NEMO code) *
30  !                                                               *
31  !************************************************************************
32  !       
33  SUBROUTINE get_partial_steps(Grid)
34    !
35    IMPLICIT NONE
36    !       
37    TYPE(Coordinates) :: Grid                     
38    REAL*8 :: za2,za1,za0,zsur,zacr,zkth,zacr2,zkth2,zdepth,zdepwp,zmin,zmax,zdiff,ze3tp,ze3wp
39    INTEGER :: i,j,jk,jj,ji,jpj,jpi,ik,ii,ipt,jpt
40    INTEGER, DIMENSION(1) :: k
41    INTEGER :: k1
42    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
43    REAL*8, POINTER, DIMENSION(:,:)   :: hdepw,e3tp,e3wp
44    REAL*8, POINTER, DIMENSION(:,:,:) :: gdept_ps,gdepw_ps
45    REAL*8 e3t_ps
46
47    !
48    WRITE(*,*) 'convert bathymetry from etopo for partial step z-coordinate case'
49    WRITE(*,*) 'minimum thickness of partial step   e3zps_min = ', e3zps_min, ' (m)'
50    WRITE(*,*) '       step  level                  e3zps_rat = ', e3zps_rat       
51    !       
52    jpi = SIZE(Grid%bathy_meter,1)
53    jpj = SIZE(Grid%bathy_meter,2)       
54    !       
55    ALLOCATE(gdepw(N),gdept(N),e3w(N),e3t(N))
56    ALLOCATE(gdepw_ps(jpi,jpj,N))                 
57    IF (.NOT.ASSOCIATED(Grid%bathy_level)) ALLOCATE(Grid%bathy_level(jpi,jpj))
58    !       
59    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
60         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
61       !   
62       WRITE(*,*) 'psur,pa0,pa1 computed'
63       za1=( ppdzmin - pphmax / (N-1) )          &
64            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
65            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
66            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
67
68       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
69       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
70       !
71    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
72         pa0.NE.0 .AND. pa1.NE.0 ) THEN
73       !       
74       WRITE(*,*) 'psur,pa0,pa1 given by namelist'
75       zsur = psur
76       za0  = pa0
77       za1  = pa1
78       za2  = pa2
79       !
80    ELSE
81       !       
82       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
83       WRITE(*,*) ' '
84       WRITE(*,*) 'please check values of variables'
85       WRITE(*,*) 'in namelist vertical_grid section'
86       WRITE(*,*) ' ' 
87       STOP   
88       !       
89    ENDIF
90
91    zacr  = ppacr
92    zkth  = ppkth     
93    zacr2 = ppacr2
94    zkth2 = ppkth2 
95    !
96    IF( ppkth == 0. ) THEN            !  uniform vertical grid
97         za1 = pphmax / FLOAT(N-1) 
98         DO i = 1, N
99            gdepw(i) = ( i - 1   ) * za1
100            gdept(i) = ( i - 0.5 ) * za1
101            e3w  (i) =  za1
102            e3t  (i) =  za1
103         END DO
104    ELSE                            ! Madec & Imbard 1996 function
105       IF( .NOT. ldbletanh ) THEN
106          DO i = 1,N
107             !
108             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
109             gdept(i) = (zsur+za0*(i+0.5)+za1*zacr*LOG(COSH(((i+0.5)-zkth)/zacr)))
110             e3w(i)   = (za0 + za1 * TANH((i-zkth)/zacr))
111             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
112             !
113          END DO
114       ELSE
115            DO i = 1,N
116               ! Double tanh function
117               gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
118                  &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
119               gdept(i) = ( zsur + za0*(i+0.5) + za1 * zacr * LOG ( COSH( ((i+0.5)-zkth ) / zacr  ) )    &
120                  &                            + za2 * zacr2* LOG ( COSH( ((i+0.5)-zkth2) / zacr2 ) )  )
121               e3w  (i) =          za0         + za1        * TANH(       (i-zkth ) / zacr  )            &
122                  &                            + za2        * TANH(       (i-zkth2) / zacr2 )
123               e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
124                  &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
125            END DO
126       ENDIF
127    ENDIF
128    gdepw(1) = 0.0 
129    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
130       !                     
131       DO i = 1, N-1
132          e3t(i) = gdepw(i+1)-gdepw(i)
133       END DO
134       e3t(N) = e3t(N-1)
135   
136       DO i = 2, N
137          e3w(i) = gdept(i) - gdept(i-1)
138       END DO
139       e3w(1  ) = 2. * (gdept(1) - gdepw(1))
140    END IF
141    !
142    ! Initialization of constant
143    !
144    zmax = gdepw(N) + e3t(N)
145    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
146    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
147    ENDIF
148    zmin = gdepw(i+1)
149    !
150    ! Initialize bathy_level to the maximum ocean level available
151    !
152    Grid%bathy_level = N-1
153    !
154    ! storage of land and island's number (zera and negative values) in mbathy
155    !
156    DO jj = 1, jpj
157       DO ji= 1, jpi
158          IF( Grid%bathy_meter(ji,jj) <= 0. )   &
159               Grid%bathy_level(ji,jj) = INT( Grid%bathy_meter(ji,jj) )
160       END DO
161    END DO
162    !
163    ! the last ocean level thickness cannot exceed e3t(jpkm1)+e3t(jpk)
164    !
165    DO jj = 1, jpj
166       DO ji= 1, jpi
167          IF( Grid%bathy_meter(ji,jj) <= 0. ) THEN
168             Grid%bathy_meter(ji,jj) = 0.e0
169          ELSE
170             Grid%bathy_meter(ji,jj) = MAX( Grid%bathy_meter(ji,jj), zmin )
171             Grid%bathy_meter(ji,jj) = MIN( Grid%bathy_meter(ji,jj), zmax )
172          ENDIF
173       END DO
174    END DO
175    !
176    ! Compute bathy_level for ocean points (i.e. the number of ocean levels)
177    ! find the number of ocean levels such that the last level thickness
178    ! is larger than the minimum of e3zps_min and e3zps_rat * e3t (where
179    ! e3t is the reference level thickness   
180    !     
181    DO jk = N-1, 1, -1
182       zdepth = gdepw(jk) + MIN( e3zps_min, e3t(jk)*e3zps_rat )
183       DO jj = 1, jpj
184          DO ji = 1, jpi
185             IF( 0. < Grid%bathy_meter(ji,jj) .AND. Grid%bathy_meter(ji,jj) <= zdepth ) &
186                  Grid%bathy_level(ji,jj) = jk-1
187          END DO
188       END DO
189    END DO
190
191
192    CALL bathymetry_control(grid%bathy_level)
193    !
194    ! initialization to the reference z-coordinate
195    !
196    WRITE(*,*) ' initialization to the reference z-coordinate '
197    !     
198    DO jk = 1, N
199       !        Write(*,*) 'k = ',jk
200       gdepw_ps(1:jpi,1:jpj,jk) = gdepw(jk)
201    END DO
202    !     
203    Grid%gdepw_ps(:,:) = gdepw_ps(:,:,3)           
204    !
205    DO jj = 1, jpj
206       DO ji = 1, jpi
207          ik = Grid%bathy_level(ji,jj)
208          ! ocean point only
209          IF( ik > 0 ) THEN
210             ! max ocean level case
211             IF( ik == N-1 ) THEN
212                zdepwp = Grid%bathy_meter(ji,jj)
213                ze3tp  = Grid%bathy_meter(ji,jj) - gdepw(ik)
214                ze3wp = 0.5 * e3w(ik) * ( 1. + ( ze3tp/e3t(ik) ) )
215                gdepw_ps(ji,jj,ik+1) = zdepwp
216                ! standard case
217             ELSE
218                !
219                IF( Grid%bathy_meter(ji,jj) <= gdepw(ik+1) ) THEN
220                   gdepw_ps(ji,jj,ik+1) = Grid%bathy_meter(ji,jj)
221                ELSE
222                   !
223                   gdepw_ps(ji,jj,ik+1) = gdepw(ik+1)
224                ENDIF
225                !
226             ENDIF
227             !           
228          ENDIF
229       END DO
230    END DO
231    !
232    DO jj = 1, jpj
233       DO ji = 1, jpi
234          ik = Grid%bathy_level(ji,jj)
235          ! ocean point only
236          IF( ik > 0 ) THEN
237             ! bathymetry output
238             !
239             Grid%gdepw_ps(ji,jj) = gdepw_ps(ji,jj,ik+1)
240             !
241             !AJOUT-----------------------------------------------------------------------
242             !
243          ELSE
244             !           
245             Grid%gdepw_ps(ji,jj) = 0
246             !
247             !AJOUT------------------------------------------------------------------------
248             !           
249          ENDIF
250          !                     
251       END DO
252    END DO
253    !     
254    !
255    DEALLOCATE(gdepw,gdept,e3w,e3t)
256    DEALLOCATE(gdepw_ps)                 
257  END SUBROUTINE get_partial_steps
258  !
259  !
260  !*************************************************************************
261  !                           *
262  ! Subroutine check interp                  *
263  !                           *
264  ! subroutine to compute gdepw_ps on the input grid (based on NEMO code) *
265  !                                                               *
266  !************************************************************************
267  !
268  !
269  SUBROUTINE check_interp( ParentGrid , gdepwChild )
270    !
271    IMPLICIT NONE
272    !                   
273    TYPE(Coordinates) :: ParentGrid
274    REAL*8,DIMENSION(:,:) :: gdepwChild 
275    INTEGER :: i,j,ji,ij,ii,jj,jpt,ipt
276    REAL,DIMENSION(N) :: gdepw,e3t
277    REAL :: za0,za1,za2,zsur,zacr,zacr2,zkth,zkth2,zmin,zmax,zdepth
278    INTEGER :: kbathy,jk,diff
279    INTEGER :: bornex,borney,bornex2,borney2
280    !
281    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
282         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
283       !   
284       WRITE(*,*) 'psur,pa0,pa1 computed'
285       za1=( ppdzmin - pphmax / (N-1) )          &
286            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
287            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
288            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
289
290       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
291       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
292       !
293    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
294         pa0.NE.0 .AND. pa1.NE.0 ) THEN
295       !       
296       WRITE(*,*) 'psur,pa0,pa1 given by namelist'
297       zsur = psur
298       za0  = pa0
299       za1  = pa1
300       za2  = pa2
301       !
302    ELSE
303       !       
304       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
305       WRITE(*,*) ' '
306       WRITE(*,*) 'please check values of variables'
307       WRITE(*,*) 'in namelist vertical_grid section'
308       WRITE(*,*) ' ' 
309       STOP   
310       !       
311    ENDIF
312
313    zacr  = ppacr
314    zkth  = ppkth     
315    zacr2 = ppacr2
316    zkth2 = ppkth2 
317    !
318    IF( ppkth == 0. ) THEN            !  uniform vertical grid
319         za1 = pphmax / FLOAT(N-1) 
320         DO i = 1, N
321            gdepw(i) = ( i - 1   ) * za1
322            e3t  (i) =  za1
323         END DO
324    ELSE                            ! Madec & Imbard 1996 function
325       IF( .NOT. ldbletanh ) THEN
326          DO i = 1,N
327             !
328             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
329             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
330             !
331          END DO
332       ELSE
333            DO i = 1,N
334               ! Double tanh function
335               gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
336                  &                       + za2 * zacr2* LOG ( COSH( (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    gdepw(1) = 0.0
343    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
344       !                     
345       DO i = 1, N-1
346          e3t(i) = gdepw(i+1)-gdepw(i)
347       END DO
348       e3t(N) = e3t(N-1)
349    END IF
350    !
351    diff = 0     
352    IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
353    !       
354    bornex = nbghostcellsfine + 1 + CEILING(irafx/2.0) + diff - irafx
355    borney = nbghostcellsfine + 1 + CEILING(irafy/2.0) + diff - irafy
356    bornex2 = nxfin - nbghostcellsfine - irafx - CEILING(irafx/2.0) 
357    borney2 = nyfin - nbghostcellsfine - irafy - CEILING(irafy/2.0)                     
358    !
359    !
360    ! west boundary
361    !
362
363    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,2+nbghostcellsfine+connectionsize*irafx-1, &
364         1,nyfin)
365
366    !
367    ! east boundary
368    !
369
370    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-1-nbghostcellsfine-(connectionsize*irafx-1),nxfin, &
371         1,nyfin)
372
373    !
374    ! north boundary
375    !
376
377    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin, &
378         nyfin-1 - nbghostcellsfine -(connectionsize*irafy-1),nyfin )
379
380    !
381    ! south boundary
382    !
383    CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin, &
384         1,2+nbghostcellsfine+connectionsize*irafy-1 )
385
386    !       
387    !
388    !
389  END SUBROUTINE check_interp
390  !
391  SUBROUTINE correct_level( gdepwchild,ParentGrid,gdepw,e3t,minboundx,maxboundx,minboundy,maxboundy )
392    !
393    IMPLICIT NONE
394    TYPE(Coordinates) :: ParentGrid
395    REAL*8,DIMENSION(:,:) :: gdepwChild
396    REAL*8,DIMENSION(N) :: gdepw,e3t
397    INTEGER :: minboundx,maxboundx,minboundy,maxboundy
398    INTEGER :: kbathy,jk,indx,indy,diff
399    REAL :: xdiff
400    INTEGER :: i,j,ji,ij,ii,jj,jpt,ipt
401    REAL*8 :: slopex, slopey,val,tmp1,tmp2,tmp3,tmp4
402    INTEGER :: parentbathy
403    REAL :: mindepth, maxdepth
404    REAL :: xmin,ymin,dxfin,dyfin,dsparent
405    INTEGER ipbegin,ipend,jpbegin,jpend
406    INTEGER ibegin,iend,jbegin,jend
407    REAL x,y,zmin,zmax
408    INTEGER ptx,pty
409    REAL,DIMENSION(:,:),ALLOCATABLE :: gdepwtemp
410    INTEGER,DIMENSION(:,:),ALLOCATABLE :: parentbathytab
411    !
412    !
413    ! Initialization of constant
414    !
415    zmax = gdepw(N) + e3t(N)
416    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
417    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
418    ENDIF
419    zmin = gdepw(i+1)
420    !
421    ! check that interpolated value stays at the same level         
422    !
423    !
424    diff = 0     
425    IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
426
427    xdiff = REAL(diff)/2.
428
429    dxfin = 1./irafx
430    dyfin = 1./irafy
431
432    ptx = 1 + nbghostcellsfine + 1
433    pty = 1 + nbghostcellsfine + 1
434
435    xmin = (imin-1) * 1
436    ymin = (jmin-1) * 1
437
438
439    ! compute x and y the locations of the indices minbounx and minboundy
440
441    x = xmin + (minboundx-ptx)*dxfin  + dxfin/2.
442    y = ymin + (minboundy-pty)*dyfin  + dyfin/2.
443
444    ! compute the indices of the nearest coarse grid points     
445    ipbegin = ptx + agrif_int((x-0.-1./2.) / 1.) - 1
446    jpbegin = pty + agrif_int((y-0.-1./2.) / 1.) - 1
447
448    ! compute indices of the fine grid points nearest to the preceeding coarse grid points       
449    ! (inferior values)
450
451    x = (ipbegin - ptx) + 1./2.
452    y = (jpbegin - pty) + 1./2.
453
454    ibegin = ptx + agrif_int((x-xmin-dxfin/2.)/dxfin) 
455    jbegin = pty + agrif_int((y-ymin-dyfin/2.)/dyfin) 
456
457    ! compute x and y the locations of the indices maxbounx and maxboundy       
458    x = xmin + (maxboundx-ptx)*dxfin + dxfin/2.
459    y = ymin + (maxboundy-pty)*dyfin + dyfin/2.
460
461    ! compute the indices of the nearest coarse grid points         
462    ipend = ptx + CEILING((x-0.-1./2) / 1.) + 1
463    jpend = pty + CEILING((y-0.-1./2) / 1.) + 1
464
465    ! compute indices of the fine grid points nearest to the preceeding coarse grid points       
466    ! (inferior values)
467
468    x = (ipend - ptx) + 1./2.
469    y = (jpend - pty) + 1./2.
470    iend = ptx + agrif_int((x-xmin-dxfin/2.)/dxfin) 
471    jend = pty + agrif_int((y-ymin-dyfin/2.)/dyfin)               
472
473
474    ALLOCATE(gdepwtemp(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
475    ALLOCATE(parentbathytab(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
476
477
478    jpt=jpbegin
479    DO j=jbegin,jend,irafy
480
481       ipt=ipbegin
482
483
484       DO i=ibegin,iend,irafx
485
486
487          !           
488          parentbathy = ParentGrid%bathy_level(ipt,jpt)
489          IF (parentbathy == 0) THEN
490             mindepth = 0.
491             maxdepth = 0.
492          ELSE
493             mindepth = MAX(gdepw(parentbathy) + MIN( e3zps_min, e3t(parentbathy)*e3zps_rat ),zmin)
494             !                  maxdepth = min(gdepw(parentbathy + 1),zmax)
495             IF (parentbathy < (N-1)) THEN
496                maxdepth = gdepw(parentbathy + 1)
497             ELSE
498                maxdepth = HUGE(1.)
499             ENDIF
500          ENDIF
501
502          slopex = vanleer(parentgrid%gdepw_ps(ipt-1:ipt+1,jpt))/REAL(irafx)
503
504
505          tmp1 = (maxdepth - parentgrid%gdepw_ps(ipt,jpt)) / REAL(irafx)
506          tmp2 = (parentgrid%gdepw_ps(ipt,jpt) - mindepth) / REAL(irafx)
507
508          IF (ABS(slopex) > tmp1) THEN
509             IF (slopex > 0) THEN
510                slopex = tmp1
511             ELSE
512                slopex = -tmp1
513             ENDIF
514          ENDIF
515
516          IF (ABS(slopex) > tmp2) THEN
517             IF (slopex > 0) THEN
518                slopex = tmp2
519             ELSE
520                slopex = -tmp2
521             ENDIF
522          ENDIF
523          !                 
524          ! interpolation on fine grid points (connection zone)
525          !
526          DO ii = i-FLOOR(irafx/2.0)+diff,i+FLOOR(irafx/2.0)
527             x = ii-i - xdiff/2.
528!!             val = parentgrid%gdepw_ps(ipt,jpt)+slopex * x
529!! chanut: uncomment this to get nearest neighbor interpolation
530             val = parentgrid%gdepw_ps(ipt,jpt)         
531             gdepwtemp(ii,j) = val
532             IF (gdepwtemp(ii,j) < mindepth) THEN
533                gdepwtemp(ii,j) = mindepth
534             ENDIF
535             IF (gdepwtemp(ii,j) > maxdepth) THEN
536                gdepwtemp(ii,j) = maxdepth
537             ENDIF
538             parentbathytab(ii,j) = parentbathy
539          ENDDO
540          ipt =ipt + 1
541       ENDDO
542
543       jpt = jpt + 1               
544    ENDDO
545
546    DO j=jbegin+irafy,jend-irafy,irafy
547
548       DO i=ibegin,iend
549
550          parentbathy = parentbathytab(i,j)
551          IF (parentbathy == 0) THEN
552             mindepth = 0.
553             maxdepth = 0.
554          ELSE
555             mindepth = MAX(gdepw(parentbathy) + MIN( e3zps_min, e3t(parentbathy)*e3zps_rat ),zmin)
556             !                  maxdepth = min(gdepw(parentbathy + 1),zmax)
557             IF (parentbathy < (N-1)) THEN
558                maxdepth = gdepw(parentbathy + 1)
559             ELSE
560                maxdepth = HUGE(1.)
561             ENDIF
562          ENDIF
563
564          slopey = vanleer(gdepwtemp(i,j-irafy:j+irafy:irafy))/REAL(irafy)
565
566          tmp1 = (maxdepth - gdepwtemp(i,j)) / REAL(irafy)
567          tmp2 = (gdepwtemp(i,j) - mindepth) / REAL(irafy)
568
569          IF (ABS(slopey) > tmp1) THEN
570             IF (slopey > 0) THEN
571                slopey = tmp1
572             ELSE
573                slopey = -tmp1
574             ENDIF
575          ENDIF
576          IF (ABS(slopey) > tmp2) THEN
577             IF (slopey > 0) THEN
578                slopey = tmp2
579             ELSE
580                slopey = -tmp2
581             ENDIF
582          ENDIF
583
584
585          DO jj = j-FLOOR(irafy/2.0)+diff,j+FLOOR(irafy/2.0)
586             y = jj-j - xdiff/2.
587!!             val = gdepwtemp(i,j) + slopey*y
588!! chanut: uncomment this to get nearest neighbor interpolation
589             val = gdepwtemp(i,j)
590             gdepwtemp(i,jj) = val     
591          ENDDO
592       ENDDO
593    ENDDO
594
595
596    gdepwchild(minboundx:maxboundx,minboundy:maxboundy) = gdepwtemp(minboundx:maxboundx,minboundy:maxboundy)
597    DEALLOCATE(gdepwtemp,parentbathytab)
598
599  END SUBROUTINE correct_level
600  !
601  !
602  !***************************************************
603  ! function van leer to compute the corresponding
604  ! Van Leer slopes
605  !***************************************************
606  !     
607  REAL FUNCTION vanleer(tab)
608    REAL, DIMENSION(3) :: tab
609    REAL res,res1
610    REAL p1,p2,p3
611
612    p1=(tab(3)-tab(1))/2.
613    p2=(tab(2)-tab(1))
614    p3=(tab(3)-tab(2))
615
616    IF ((p1>0.).AND.(p2>0.).AND.(p3>0)) THEN
617       res1=MINVAL((/p1,p2,p3/))
618    ELSEIF ((p1<0.).AND.(p2<0.).AND.(p3<0)) THEN
619       res1=MAXVAL((/p1,p2,p3/))
620    ELSE
621       res1=0.
622    ENDIF
623
624    vanleer = res1   
625
626
627  END FUNCTION vanleer
628  !
629  !
630  !********************************************************************************
631  !   subroutine bathymetry_control                               *
632  !                                                            *
633  !  - Purpose :   check the bathymetry in levels                       *
634  !                              *
635  !  - Method  :   The array mbathy is checked to verified its consistency *
636  !      with the model options. in particular:             *
637  !            mbathy must have at least 1 land grid-points (mbathy<=0)    *
638  !                  along closed boundary.              *
639  !            mbathy must be cyclic IF jperio=1.              *
640  !            mbathy must be lower or equal to jpk-1.            *
641  !            isolated ocean grid points are suppressed from mbathy    *
642  !                  since they are only connected to remaining         *
643  !                  ocean through vertical diffusion.            *
644  !                              *
645  !                              *
646  !********************************************************************************
647
648  SUBROUTINE bathymetry_control(mbathy)
649
650    INTEGER ::   i, j, jl           
651    INTEGER ::   icompt, ibtest, ikmax         
652    REAL*8, DIMENSION(:,:) :: mbathy     
653
654    ! ================
655    ! Bathymetry check
656    ! ================
657
658    ! Suppress isolated ocean grid points
659
660    WRITE(*,*)'                   suppress isolated ocean grid points'
661    WRITE(*,*)'                   -----------------------------------'
662
663    icompt = 0
664
665    DO jl = 1, 2
666       !
667       DO j = 2, SIZE(mbathy,2)-1
668          DO i = 2, SIZE(mbathy,1)-1
669
670             ibtest = MAX( mbathy(i-1,j), mbathy(i+1,j),mbathy(i,j-1),mbathy(i,j+1) )
671             !               
672             IF( ibtest < mbathy(i,j) ) THEN
673                !                 
674                WRITE(*,*) 'grid-point(i,j)= ',i,j,'is changed from',mbathy(i,j),' to ', ibtest
675                mbathy(i,j) = ibtest
676                icompt = icompt + 1
677                !
678             ENDIF
679             !           
680          END DO
681       END DO
682       !
683    END DO
684    !     
685    IF( icompt == 0 ) THEN
686       WRITE(*,*)'     no isolated ocean grid points'
687    ELSE
688       WRITE(*,*)'    ',icompt,' ocean grid points suppressed'
689    ENDIF
690    !
691
692    ! Number of ocean level inferior or equal to jpkm1
693
694    ikmax = 0
695    DO j = 1, SIZE(mbathy,2)
696       DO ji = 1, SIZE(mbathy,1)
697          ikmax = MAX( ikmax, NINT(mbathy(i,j)) )
698       END DO
699    END DO
700    !
701    IF( ikmax > N-1 ) THEN
702       WRITE(*,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
703       WRITE(*,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
704    ELSE IF( ikmax < N-1 ) THEN
705       WRITE(*,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
706       WRITE(*,*) ' you can decrease jpk to ', ikmax+1
707    ENDIF
708
709  END SUBROUTINE bathymetry_control
710  !
711  !
712  !**********************************************************************************
713  !
714  !subroutine get_scale_factors
715  !
716  !**********************************************************************************
717  !
718  SUBROUTINE get_scale_factors(Grid,fse3t,fse3u,fse3v)
719    !
720    IMPLICIT NONE
721    !       
722    TYPE(Coordinates) :: Grid 
723    REAL*8, DIMENSION(:,:,:) :: fse3u,fse3t,fse3v
724    !                                 
725    REAL*8 :: za2,za1,za0,zsur,zacr,zkth,zacr2,zkth2,zdepth,zdepwp,zmin,zmax,zdiff,ze3tp,ze3wp
726    INTEGER :: i,j,jk,jj,ji,jpj,jpi,ik,ii,ipt,jpt,jpk
727    INTEGER, DIMENSION(1) :: k
728    INTEGER :: k1
729    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
730    REAL*8, POINTER, DIMENSION(:,:)   :: hdepw,e3tp,e3wp
731    REAL*8, POINTER, DIMENSION(:,:,:) :: gdept_ps,gdepw_ps   
732    !       
733    jpi = SIZE(fse3t,1)
734    jpj = SIZE(fse3t,2) 
735    jpk = SIZE(fse3t,3)     
736    !       
737    ALLOCATE(gdepw(jpk),e3t(jpk))
738    ALLOCATE(gdepw_ps(jpi,jpj,jpk))                 
739   
740    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
741         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
742       !   
743       WRITE(*,*) 'psur,pa0,pa1 computed'
744       za1=( ppdzmin - pphmax / (jpk-1) )          &
745            / ( TANH((1-ppkth)/ppacr) - ppacr/(jpk-1) &
746            *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
747            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
748
749       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
750       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
751       !
752    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
753         pa0.NE.0 .AND. pa1.NE.0 ) THEN
754       !       
755       WRITE(*,*) 'psur,pa0,pa1 given by namelist'
756       zsur = psur
757       za0  = pa0
758       za1  = pa1
759       za2  = pa2
760       !
761    ELSE
762       !       
763       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
764       WRITE(*,*) ' '
765       WRITE(*,*) 'please check values of variables'
766       WRITE(*,*) 'in namelist vertical_grid section'
767       WRITE(*,*) ' ' 
768       STOP   
769       !       
770    ENDIF
771
772    zacr  = ppacr
773    zkth  = ppkth     
774    zacr2 = ppacr2
775    zkth2 = ppkth2 
776    !
777    IF( ppkth == 0. ) THEN            !  uniform vertical grid
778         za1 = pphmax / FLOAT(jpk-1) 
779         DO i = 1, jpk
780            gdepw(i) = ( i - 1   ) * za1
781            e3t  (i) =  za1
782         END DO
783    ELSE                            ! Madec & Imbard 1996 function
784       IF( .NOT. ldbletanh ) THEN
785          DO i = 1,jpk
786             !
787             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
788             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
789             !
790          END DO
791       ELSE
792            DO i = 1,jpk
793               ! Double tanh function
794               gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
795                  &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
796               e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
797                  &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
798            END DO
799       ENDIF
800    ENDIF         
801    !         
802    gdepw(1)=0.
803    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
804       !                     
805       DO i = 1, jpk-1
806          e3t(i) = gdepw(i+1)-gdepw(i)
807       END DO
808       e3t(jpk) = e3t(jpk-1)
809    END IF
810    !               
811    DO i = 1,jpk
812       !       
813       fse3t(:,:,i) = e3t(i)
814       gdepw_ps(:,:,i) = gdepw(i)
815       !
816    END DO
817    !                                 
818    gdepw(1) = 0.0 
819    gdepw_ps(:,:,1) = 0.0 
820    !
821    zmax = gdepw(jpk) + e3t(jpk)
822    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
823    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
824    ENDIF
825    zmin = gdepw(i+1)
826    !
827    DO jj = 1, jpj
828       DO ji= 1, jpi
829          IF( Grid%bathy_meter(ji,jj) <= 0. ) THEN
830             Grid%bathy_meter(ji,jj) = 0.e0
831          ELSE
832             Grid%bathy_meter(ji,jj) = MAX( Grid%bathy_meter(ji,jj), zmin )
833             Grid%bathy_meter(ji,jj) = MIN( Grid%bathy_meter(ji,jj), zmax )
834          ENDIF
835       END DO
836    END DO
837    !
838    DO jj = 1, jpj
839       DO ji = 1, jpi
840          ik = Grid%bathy_level(ji,jj) 
841          IF( ik > 0 ) THEN
842             ! max ocean level case
843             IF( ik == jpk-1 ) THEN
844                zdepwp = Grid%bathy_meter(ji,jj)
845                ze3tp  = Grid%bathy_meter(ji,jj) - gdepw(ik)
846                fse3t(ji,jj,ik  ) = ze3tp
847                fse3t(ji,jj,ik+1) = ze3tp
848                gdepw_ps(ji,jj,ik+1) = zdepwp
849             ELSE
850                IF( Grid%bathy_meter(ji,jj) <= gdepw(ik+1) ) THEN
851                   gdepw_ps(ji,jj,ik+1) = Grid%bathy_meter(ji,jj)
852                ELSE
853                   gdepw_ps(ji,jj,ik+1) = gdepw(ik+1)
854                ENDIF
855                fse3t(ji,jj,ik) = e3t(ik) * ( gdepw_ps(ji,jj,ik+1) - gdepw(ik))   & 
856                     /( gdepw(ik+1) - gdepw(ik)) 
857                fse3t(ji,jj,ik+1) = fse3t(ji,jj,ik)
858
859             ENDIF
860          ENDIF
861       END DO
862    END DO
863    !
864    DO i = 1, jpk
865       fse3u (:,:,i)  = e3t(i)
866       fse3v (:,:,i)  = e3t(i)
867    END DO
868    !
869    DO jk = 1,jpk
870       DO jj = 1, jpj-1
871          DO ji = 1, jpi-1
872             fse3u (ji,jj,jk) = MIN( fse3t(ji,jj,jk), fse3t(ji+1,jj,jk))
873             fse3v (ji,jj,jk) = MIN( fse3t(ji,jj,jk), fse3t(ji,jj+1,jk))     
874          ENDDO
875       ENDDO
876    ENDDO
877    !           
878    DEALLOCATE(gdepw,e3t)
879    DEALLOCATE(gdepw_ps)
880    DEALLOCATE(Grid%bathy_meter,Grid%bathy_level) 
881    !
882  END SUBROUTINE get_scale_factors
883  !       
884END MODULE agrif_partial_steps
885
886
Note: See TracBrowser for help on using the repository browser.