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

source: branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/TOOLS/NESTING/src/agrif_partial_steps.f90 @ 5641

Last change on this file since 5641 was 5641, checked in by jchanut, 9 years ago

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