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

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

Last change on this file since 5656 was 5656, checked in by timgraham, 9 years ago

Merge of AGRIF branch (branches/2014/dev_r4765_CNRS_agrif) onto the trunk
  • 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    zmin = gdepw(i+1)
401    !
402    ! check that interpolated value stays at the same level         
403    !
404    !
405    diff = 0     
406    IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
407
408    xdiff = REAL(diff)/2.
409
410    dxfin = 1./irafx
411    dyfin = 1./irafy
412
413    ptx = 3
414    pty = 3
415
416    xmin = (imin-1) * 1
417    ymin = (jmin-1) * 1
418
419
420    ! compute x and y the locations of the indices minbounx and minboundy
421
422    x = xmin + (minboundx-ptx)*dxfin  + dxfin/2.
423    y = ymin + (minboundy-pty)*dyfin  + dyfin/2.
424
425    ! compute the indices of the nearest coarse grid points     
426    ipbegin = ptx + agrif_int((x-0.-1./2.) / 1.) - 1
427    jpbegin = pty + agrif_int((y-0.-1./2.) / 1.) - 1
428
429    ! compute indices of the fine grid points nearest to the preceeding coarse grid points       
430    ! (inferior values)
431
432    x = (ipbegin - ptx) + 1./2.
433    y = (jpbegin - pty) + 1./2.
434
435    ibegin = ptx + agrif_int((x-xmin-dxfin/2.)/dxfin) 
436    jbegin = pty + agrif_int((y-ymin-dyfin/2.)/dyfin) 
437
438    ! compute x and y the locations of the indices maxbounx and maxboundy       
439    x = xmin + (maxboundx-ptx)*dxfin + dxfin/2.
440    y = ymin + (maxboundy-pty)*dyfin + dyfin/2.
441
442    ! compute the indices of the nearest coarse grid points         
443    ipend = ptx + CEILING((x-0.-1./2) / 1.) + 1
444    jpend = pty + CEILING((y-0.-1./2) / 1.) + 1
445
446    ! compute indices of the fine grid points nearest to the preceeding coarse grid points       
447    ! (inferior values)
448
449    x = (ipend - ptx) + 1./2.
450    y = (jpend - pty) + 1./2.
451    iend = ptx + agrif_int((x-xmin-dxfin/2.)/dxfin) 
452    jend = pty + agrif_int((y-ymin-dyfin/2.)/dyfin)               
453
454
455    ALLOCATE(gdepwtemp(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
456    ALLOCATE(parentbathytab(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
457
458
459    jpt=jpbegin
460    DO j=jbegin,jend,irafy
461
462       ipt=ipbegin
463
464
465       DO i=ibegin,iend,irafx
466
467
468          !           
469          parentbathy = ParentGrid%bathy_level(ipt,jpt)
470          IF (parentbathy == 0) THEN
471             mindepth = 0.
472             maxdepth = 0.
473          ELSE
474             mindepth = MAX(gdepw(parentbathy) + MIN( e3zps_min, e3t(parentbathy)*e3zps_rat ),zmin)
475             !                  maxdepth = min(gdepw(parentbathy + 1),zmax)
476             IF (parentbathy < (N-1)) THEN
477                maxdepth = gdepw(parentbathy + 1)
478             ELSE
479                maxdepth = HUGE(1.)
480             ENDIF
481          ENDIF
482
483          slopex = vanleer(parentgrid%gdepw_ps(ipt-1:ipt+1,jpt))/REAL(irafx)
484
485
486          tmp1 = (maxdepth - parentgrid%gdepw_ps(ipt,jpt)) / REAL(irafx)
487          tmp2 = (parentgrid%gdepw_ps(ipt,jpt) - mindepth) / REAL(irafx)
488
489          IF (ABS(slopex) > tmp1) THEN
490             IF (slopex > 0) THEN
491                slopex = tmp1
492             ELSE
493                slopex = -tmp1
494             ENDIF
495          ENDIF
496
497          IF (ABS(slopex) > tmp2) THEN
498             IF (slopex > 0) THEN
499                slopex = tmp2
500             ELSE
501                slopex = -tmp2
502             ENDIF
503          ENDIF
504          !                 
505          ! interpolation on fine grid points (connection zone)
506          !
507          DO ii = i-FLOOR(irafx/2.0)+diff,i+FLOOR(irafx/2.0)
508             x = ii-i - xdiff/2.
509             val = parentgrid%gdepw_ps(ipt,jpt)+slopex * x
510             gdepwtemp(ii,j) = val
511             IF (gdepwtemp(ii,j) < mindepth) THEN
512                gdepwtemp(ii,j) = mindepth
513             ENDIF
514             IF (gdepwtemp(ii,j) > maxdepth) THEN
515                gdepwtemp(ii,j) = maxdepth
516             ENDIF
517             parentbathytab(ii,j) = parentbathy
518          ENDDO
519          ipt =ipt + 1
520       ENDDO
521
522       jpt = jpt + 1               
523    ENDDO
524
525    DO j=jbegin+irafy,jend-irafy,irafy
526
527       DO i=ibegin,iend
528
529          parentbathy = parentbathytab(i,j)
530          IF (parentbathy == 0) THEN
531             mindepth = 0.
532             maxdepth = 0.
533          ELSE
534             mindepth = MAX(gdepw(parentbathy) + MIN( e3zps_min, e3t(parentbathy)*e3zps_rat ),zmin)
535             !                  maxdepth = min(gdepw(parentbathy + 1),zmax)
536             IF (parentbathy < (N-1)) THEN
537                maxdepth = gdepw(parentbathy + 1)
538             ELSE
539                maxdepth = HUGE(1.)
540             ENDIF
541          ENDIF
542
543          slopey = vanleer(gdepwtemp(i,j-irafy:j+irafy:irafy))/REAL(irafy)
544
545          tmp1 = (maxdepth - gdepwtemp(i,j)) / REAL(irafy)
546          tmp2 = (gdepwtemp(i,j) - mindepth) / REAL(irafy)
547
548          IF (ABS(slopey) > tmp1) THEN
549             IF (slopey > 0) THEN
550                slopey = tmp1
551             ELSE
552                slopey = -tmp1
553             ENDIF
554          ENDIF
555          IF (ABS(slopey) > tmp2) THEN
556             IF (slopey > 0) THEN
557                slopey = tmp2
558             ELSE
559                slopey = -tmp2
560             ENDIF
561          ENDIF
562
563
564          DO jj = j-FLOOR(irafy/2.0)+diff,j+FLOOR(irafy/2.0)
565             y = jj-j - xdiff/2.
566             val = gdepwtemp(i,j) + slopey*y
567             gdepwtemp(i,jj) = val     
568          ENDDO
569       ENDDO
570    ENDDO
571
572
573    gdepwchild(minboundx:maxboundx,minboundy:maxboundy) = gdepwtemp(minboundx:maxboundx,minboundy:maxboundy)
574    DEALLOCATE(gdepwtemp,parentbathytab)
575
576  END SUBROUTINE correct_level
577  !
578  !
579  !***************************************************
580  ! function van leer to compute the corresponding
581  ! Van Leer slopes
582  !***************************************************
583  !     
584  REAL FUNCTION vanleer(tab)
585    REAL, DIMENSION(3) :: tab
586    REAL res,res1
587    REAL p1,p2,p3
588
589    p1=(tab(3)-tab(1))/2.
590    p2=(tab(2)-tab(1))
591    p3=(tab(3)-tab(2))
592
593    IF ((p1>0.).AND.(p2>0.).AND.(p3>0)) THEN
594       res1=MINVAL((/p1,p2,p3/))
595    ELSEIF ((p1<0.).AND.(p2<0.).AND.(p3<0)) THEN
596       res1=MAXVAL((/p1,p2,p3/))
597    ELSE
598       res1=0.
599    ENDIF
600
601    vanleer = res1   
602
603
604  END FUNCTION vanleer
605  !
606  !
607  !********************************************************************************
608  !   subroutine bathymetry_control                               *
609  !                                                            *
610  !  - Purpose :   check the bathymetry in levels                       *
611  !                              *
612  !  - Method  :   The array mbathy is checked to verified its consistency *
613  !      with the model options. in particular:             *
614  !            mbathy must have at least 1 land grid-points (mbathy<=0)    *
615  !                  along closed boundary.              *
616  !            mbathy must be cyclic IF jperio=1.              *
617  !            mbathy must be lower or equal to jpk-1.            *
618  !            isolated ocean grid points are suppressed from mbathy    *
619  !                  since they are only connected to remaining         *
620  !                  ocean through vertical diffusion.            *
621  !                              *
622  !                              *
623  !********************************************************************************
624
625  SUBROUTINE bathymetry_control(mbathy)
626
627    INTEGER ::   i, j, jl           
628    INTEGER ::   icompt, ibtest, ikmax         
629    REAL*8, DIMENSION(:,:) :: mbathy     
630
631    ! ================
632    ! Bathymetry check
633    ! ================
634
635    ! Suppress isolated ocean grid points
636
637    WRITE(*,*)'                   suppress isolated ocean grid points'
638    WRITE(*,*)'                   -----------------------------------'
639
640    icompt = 0
641
642    DO jl = 1, 2
643       !
644       DO j = 2, SIZE(mbathy,2)-1
645          DO i = 2, SIZE(mbathy,1)-1
646
647             ibtest = MAX( mbathy(i-1,j), mbathy(i+1,j),mbathy(i,j-1),mbathy(i,j+1) )
648             !               
649             IF( ibtest < mbathy(i,j) ) THEN
650                !                 
651                WRITE(*,*) 'grid-point(i,j)= ',i,j,'is changed from',mbathy(i,j),' to ', ibtest
652                mbathy(i,j) = ibtest
653                icompt = icompt + 1
654                !
655             ENDIF
656             !           
657          END DO
658       END DO
659       !
660    END DO
661    !     
662    IF( icompt == 0 ) THEN
663       WRITE(*,*)'     no isolated ocean grid points'
664    ELSE
665       WRITE(*,*)'    ',icompt,' ocean grid points suppressed'
666    ENDIF
667    !
668
669    ! Number of ocean level inferior or equal to jpkm1
670
671    ikmax = 0
672    DO j = 1, SIZE(mbathy,2)
673       DO ji = 1, SIZE(mbathy,1)
674          ikmax = MAX( ikmax, NINT(mbathy(i,j)) )
675       END DO
676    END DO
677    !
678    IF( ikmax > N-1 ) THEN
679       WRITE(*,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
680       WRITE(*,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
681    ELSE IF( ikmax < N-1 ) THEN
682       WRITE(*,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
683       WRITE(*,*) ' you can decrease jpk to ', ikmax+1
684    ENDIF
685
686  END SUBROUTINE bathymetry_control
687  !
688  !
689  !**********************************************************************************
690  !
691  !subroutine get_scale_factors
692  !
693  !**********************************************************************************
694  !
695  SUBROUTINE get_scale_factors(Grid,fse3t,fse3u,fse3v)
696    !
697    IMPLICIT NONE
698    !       
699    TYPE(Coordinates) :: Grid 
700    REAL*8, DIMENSION(:,:,:) :: fse3u,fse3t,fse3v
701    !                                 
702    REAL*8 :: za2,za1,za0,zsur,zacr,zkth,zacr2,zkth2,zdepth,zdepwp,zmin,zmax,zdiff,ze3tp,ze3wp
703    INTEGER :: i,j,jk,jj,ji,jpj,jpi,ik,ii,ipt,jpt,jpk
704    INTEGER, DIMENSION(1) :: k
705    INTEGER :: k1
706    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
707    REAL*8, POINTER, DIMENSION(:,:)   :: hdepw,e3tp,e3wp
708    REAL*8, POINTER, DIMENSION(:,:,:) :: gdept_ps,gdepw_ps   
709    !       
710    jpi = SIZE(fse3t,1)
711    jpj = SIZE(fse3t,2) 
712    jpk = SIZE(fse3t,3)     
713    !       
714    ALLOCATE(gdepw(jpk),e3t(jpk))
715    ALLOCATE(gdepw_ps(jpi,jpj,jpk))                 
716    ! 
717    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
718         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
719       !   
720       WRITE(*,*) 'psur,pa0,pa1 computed'
721       za1=( ppdzmin - pphmax / (jpk-1) )          &
722            / ( TANH((1-ppkth)/ppacr) - ppacr/(jpk-1) &
723            *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
724            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
725
726       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
727       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
728       !
729    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
730         pa0.NE.0 .AND. pa1.NE.0 ) THEN
731       !       
732       WRITE(*,*) 'psur,pa0,pa1 given by namelist'
733       zsur = psur
734       za0  = pa0
735       za1  = pa1
736       za2  = pa2
737       !
738    ELSE
739       !       
740       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
741       WRITE(*,*) ' '
742       WRITE(*,*) 'please check values of variables'
743       WRITE(*,*) 'in namelist vertical_grid section'
744       WRITE(*,*) ' ' 
745       STOP   
746       !       
747    ENDIF
748
749    zacr  = ppacr
750    zkth  = ppkth     
751    zacr2 = ppacr2
752    zkth2 = ppkth2 
753    !
754    IF( ppkth == 0. ) THEN            !  uniform vertical grid
755         za1 = pphmax / FLOAT(jpk-1) 
756         DO i = 1, jpk
757            gdepw(i) = ( i - 1   ) * za1
758            e3t  (i) =  za1
759         END DO
760    ELSE                            ! Madec & Imbard 1996 function
761       IF( .NOT. ldbletanh ) THEN
762          DO i = 1,jpk
763             !
764             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
765             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
766             !
767          END DO
768       ELSE
769            DO i = 1,jpk
770               ! Double tanh function
771               gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
772                  &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
773               e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
774                  &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
775            END DO
776       ENDIF
777    ENDIF         
778    !         
779    !               
780    DO i = 1,jpk
781       !       
782       fse3t(:,:,i) = e3t(i)
783       gdepw_ps(:,:,i) = gdepw(i)
784       !
785    END DO
786    !                                 
787    gdepw(1) = 0.0 
788    gdepw_ps(:,:,1) = 0.0 
789    !
790    zmax = gdepw(jpk) + e3t(jpk)
791    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
792    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
793    ENDIF
794    zmin = gdepw(i+1)
795    !
796    DO jj = 1, jpj
797       DO ji= 1, jpi
798          IF( Grid%bathy_meter(ji,jj) <= 0. ) THEN
799             Grid%bathy_meter(ji,jj) = 0.e0
800          ELSE
801             Grid%bathy_meter(ji,jj) = MAX( Grid%bathy_meter(ji,jj), zmin )
802             Grid%bathy_meter(ji,jj) = MIN( Grid%bathy_meter(ji,jj), zmax )
803          ENDIF
804       END DO
805    END DO
806    !
807    DO jj = 1, jpj
808       DO ji = 1, jpi
809          ik = Grid%bathy_level(ji,jj) 
810          IF( ik > 0 ) THEN
811             ! max ocean level case
812             IF( ik == jpk-1 ) THEN
813                zdepwp = Grid%bathy_meter(ji,jj)
814                ze3tp  = Grid%bathy_meter(ji,jj) - gdepw(ik)
815                fse3t(ji,jj,ik  ) = ze3tp
816                fse3t(ji,jj,ik+1) = ze3tp
817                gdepw_ps(ji,jj,ik+1) = zdepwp
818             ELSE
819                IF( Grid%bathy_meter(ji,jj) <= gdepw(ik+1) ) THEN
820                   gdepw_ps(ji,jj,ik+1) = Grid%bathy_meter(ji,jj)
821                ELSE
822                   gdepw_ps(ji,jj,ik+1) = gdepw(ik+1)
823                ENDIF
824                fse3t(ji,jj,ik) = e3t(ik) * ( gdepw_ps(ji,jj,ik+1) - gdepw(ik))   & 
825                     /( gdepw(ik+1) - gdepw(ik)) 
826                fse3t(ji,jj,ik+1) = fse3t(ji,jj,ik)
827
828             ENDIF
829          ENDIF
830       END DO
831    END DO
832    !
833    DO i = 1, jpk
834       fse3u (:,:,i)  = e3t(i)
835       fse3v (:,:,i)  = e3t(i)
836    END DO
837    !
838    DO jk = 1,jpk
839       DO jj = 1, jpj-1
840          DO ji = 1, jpi-1
841             fse3u (ji,jj,jk) = MIN( fse3t(ji,jj,jk), fse3t(ji+1,jj,jk))
842             fse3v (ji,jj,jk) = MIN( fse3t(ji,jj,jk), fse3t(ji,jj+1,jk))     
843          ENDDO
844       ENDDO
845    ENDDO
846    !           
847    DEALLOCATE(gdepw,e3t)
848    DEALLOCATE(gdepw_ps)
849    DEALLOCATE(Grid%bathy_meter,Grid%bathy_level) 
850    !
851  END SUBROUTINE get_scale_factors
852  !       
853END MODULE agrif_partial_steps
854
855
Note: See TracBrowser for help on using the repository browser.