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 utils/tools/NESTING/src – NEMO

source: utils/tools/NESTING/src/agrif_partial_steps.f90 @ 10025

Last change on this file since 10025 was 10025, checked in by clem, 6 years ago

nesting tools are partly rewritten (mostly for create_coordinates and bathy) to get better functionality. Now you can use the nesting to either define an agrif zoom or a regional domain (for bdy purposes). Also, the nesting tools output a domain_cfg.nc that can be directly used in NEMO4 without the need of DOMAINcfg tool. The option of median average for bathymetry interpolation still does not work properly but it's not new
  • Property svn:keywords set to Id
File size: 29.5 KB
Line 
1!************************************************************************
2! Fortran 95 OPA Nesting tools                                       *
3!                                                              *
4!     Copyright (C) 2005 Florian Lemarié (Florian.Lemarie@imag.fr)      *
5!                        Laurent Debreu (Laurent.Debreu@imag.fr)        *
6!************************************************************************
7!
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
279    !
280    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
281         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
282       !   
283       WRITE(*,*) 'psur,pa0,pa1 computed'
284       za1=( ppdzmin - pphmax / (N-1) )          &
285            / ( TANH((1-ppkth)/ppacr) - ppacr/(N-1) &
286            *  (  LOG( COSH( (N - ppkth) / ppacr) )      &
287            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
288
289       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
290       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
291       !
292    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
293         pa0.NE.0 .AND. pa1.NE.0 ) THEN
294       !       
295       WRITE(*,*) 'psur,pa0,pa1 given by namelist'
296       zsur = psur
297       za0  = pa0
298       za1  = pa1
299       za2  = pa2
300       !
301    ELSE
302       !       
303       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
304       WRITE(*,*) ' '
305       WRITE(*,*) 'please check values of variables'
306       WRITE(*,*) 'in namelist vertical_grid section'
307       WRITE(*,*) ' ' 
308       STOP   
309       !       
310    ENDIF
311
312    zacr  = ppacr
313    zkth  = ppkth     
314    zacr2 = ppacr2
315    zkth2 = ppkth2 
316    !
317    IF( ppkth == 0. ) THEN            !  uniform vertical grid
318         za1 = pphmax / FLOAT(N-1) 
319         DO i = 1, N
320            gdepw(i) = ( i - 1   ) * za1
321            e3t  (i) =  za1
322         END DO
323    ELSE                            ! Madec & Imbard 1996 function
324       IF( .NOT. ldbletanh ) THEN
325          DO i = 1,N
326             !
327             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
328             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
329             !
330          END DO
331       ELSE
332            DO i = 1,N
333               ! Double tanh function
334               gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
335                  &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
336               e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
337                  &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
338            END DO
339       ENDIF
340    ENDIF
341    gdepw(1) = 0.0
342    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
343       !                     
344       DO i = 1, N-1
345          e3t(i) = gdepw(i+1)-gdepw(i)
346       END DO
347       e3t(N) = e3t(N-1)
348    END IF
349    !             
350    !
351    ! west boundary
352    IF( ln_agrif_domain ) THEN
353       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafx-1,1,nyfin)
354    ELSE
355       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,(npt_copy+npt_connect)*irafx,1,nyfin)
356    ENDIF
357    !
358    ! east boundary
359    IF( ln_agrif_domain ) THEN
360       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafx-1),nxfin,1,nyfin)
361    ELSE
362       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,nxfin-((npt_copy+npt_connect)*irafx+1),nxfin,1,nyfin)
363    ENDIF
364    !
365    ! north boundary
366    IF( ln_agrif_domain ) THEN
367       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,nyfin-1-nbghostcellsfine-((npt_copy+npt_connect)*irafy-1),nyfin)
368    ELSE
369       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,nyfin-((npt_copy+npt_connect)*irafy+1),nyfin)
370    ENDIF
371    !
372    ! south boundary
373    IF( ln_agrif_domain ) THEN
374       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,1,2+nbghostcellsfine+(npt_copy+npt_connect)*irafy-1)
375    ELSE
376       CALL correct_level( gdepwchild,ParentGrid,gdepw,e3t,1,nxfin,1,(npt_copy+npt_connect)*irafy)
377    ENDIF
378    !       
379    !
380    !
381  END SUBROUTINE check_interp
382  !
383  SUBROUTINE correct_level( gdepwchild,ParentGrid,gdepw,e3t,minboundx,maxboundx,minboundy,maxboundy )
384    !
385    IMPLICIT NONE
386    TYPE(Coordinates) :: ParentGrid
387    REAL*8,DIMENSION(:,:) :: gdepwChild
388    REAL*8,DIMENSION(N) :: gdepw,e3t
389    INTEGER :: minboundx,maxboundx,minboundy,maxboundy
390    INTEGER :: kbathy,jk,indx,indy,diff
391    REAL :: xdiff
392    INTEGER :: i,j,ji,ij,ii,jj,jpt,ipt,i1,i2,j1,j2,ii1,ii2,jj1,jj2
393    REAL*8 :: slopex, slopey,val,tmp1,tmp2,tmp3,tmp4
394    INTEGER :: parentbathy
395    REAL :: mindepth, maxdepth
396    REAL :: xmin,ymin,dxfin,dyfin,dsparent
397    INTEGER ipbegin,ipend,jpbegin,jpend
398    INTEGER ibegin,iend,jbegin,jend
399    REAL x,y,zmin,zmax
400    INTEGER ptx,pty
401    REAL,DIMENSION(:,:),ALLOCATABLE :: gdepwtemp
402    INTEGER,DIMENSION(:,:),ALLOCATABLE :: parentbathytab
403    !
404    !
405    ! Initialization of constant
406    !
407    zmax = gdepw(N) + e3t(N)
408    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
409    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
410    ENDIF
411    zmin = gdepw(i+1)
412    !
413    ! check that interpolated value stays at the same level         
414    !
415    !
416    diff = 0     
417    IF ( MOD(irafx,2) .EQ. 0 ) diff = 1
418
419    xdiff = REAL(diff)/2.
420
421    dxfin = 1./irafx
422    dyfin = 1./irafy
423
424    ptx = 1 + nbghostcellsfine + 1
425    pty = 1 + nbghostcellsfine + 1
426
427    xmin = (imin-1) * 1
428    ymin = (jmin-1) * 1
429
430
431    ! compute x and y the locations of the indices minbounx and minboundy
432
433    x = xmin + (minboundx-ptx)*dxfin  + dxfin/2.
434    y = ymin + (minboundy-pty)*dyfin  + dyfin/2.
435
436    ! compute the indices of the nearest coarse grid points     
437    ipbegin = ptx + agrif_int((x-0.-1./2.) / 1.) - 1
438    jpbegin = pty + agrif_int((y-0.-1./2.) / 1.) - 1
439
440    ! compute indices of the fine grid points nearest to the preceeding coarse grid points       
441    ! (inferior values)
442
443    x = (ipbegin - ptx) + 1./2.
444    y = (jpbegin - pty) + 1./2.
445
446    ibegin = ptx + agrif_int((x-xmin-dxfin/2.)/dxfin) 
447    jbegin = pty + agrif_int((y-ymin-dyfin/2.)/dyfin) 
448
449    ! compute x and y the locations of the indices maxbounx and maxboundy       
450    x = xmin + (maxboundx-ptx)*dxfin + dxfin/2.
451    y = ymin + (maxboundy-pty)*dyfin + dyfin/2.
452
453    ! compute the indices of the nearest coarse grid points         
454    ipend = ptx + CEILING((x-0.-1./2) / 1.) + 1
455    jpend = pty + CEILING((y-0.-1./2) / 1.) + 1
456
457    ! compute indices of the fine grid points nearest to the preceeding coarse grid points       
458    ! (inferior values)
459
460    x = (ipend - ptx) + 1./2.
461    y = (jpend - pty) + 1./2.
462    iend = ptx + agrif_int((x-xmin-dxfin/2.)/dxfin) 
463    jend = pty + agrif_int((y-ymin-dyfin/2.)/dyfin)               
464
465    IF( ln_agrif_domain ) THEN
466       ALLOCATE(gdepwtemp(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
467       ALLOCATE(parentbathytab(ibegin-irafx:iend+irafx,jbegin-irafy:jend+irafy))
468
469       i1 = ibegin
470       i2 = iend
471       j1 = jbegin
472       j2 = jend
473       
474       ii1 = -FLOOR(irafx/2.0)+diff
475       ii2 =  FLOOR(irafx/2.0)
476       jj1 = -FLOOR(irafy/2.0)+diff
477       jj2 =  FLOOR(irafy/2.0)
478    ELSE
479       ibegin = minboundx
480       jbegin = minboundy
481       iend   = maxboundx ! (npt_copy+npt_connect)*irafx
482       jend   = maxboundy ! (npt_copy+npt_connect)*irafy
483       !
484       ipbegin = imin + (ibegin-1)/irafx
485       jpbegin = jmin + (jbegin-1)/irafy
486       ipend   = ipbegin + (npt_copy+npt_connect) - 1
487       jpend   = jpbegin + (npt_copy+npt_connect) - 1
488       !
489       i1 = ibegin
490       i2 = iend
491       j1 = jbegin
492       j2 = jend
493       
494       ii1 = 0
495       ii2 = irafx - 1
496       jj1 = 0
497       jj2 = irafy - 1
498       !
499       ALLOCATE(gdepwtemp(ibegin:iend,jbegin:jend))
500       ALLOCATE(parentbathytab(ibegin:iend,jbegin:jend))
501
502    ENDIF
503   
504
505    jpt=jpbegin
506    DO j=jbegin,jend,irafy
507
508       ipt=ipbegin
509
510
511       DO i=i1,i2,irafx
512
513
514          !           
515          parentbathy = ParentGrid%bathy_level(ipt,jpt)
516          IF (parentbathy == 0) THEN
517             mindepth = 0.
518             maxdepth = 0.
519          ELSE
520             mindepth = MAX(gdepw(parentbathy) + MIN( e3zps_min, e3t(parentbathy)*e3zps_rat ),zmin)
521             !                  maxdepth = min(gdepw(parentbathy + 1),zmax)
522             IF (parentbathy < (N-1)) THEN
523                maxdepth = gdepw(parentbathy + 1)
524             ELSE
525                maxdepth = HUGE(1.)
526             ENDIF
527          ENDIF
528
529          slopex = vanleer(parentgrid%gdepw_ps(ipt-1:ipt+1,jpt))/REAL(irafx)
530
531
532          tmp1 = (maxdepth - parentgrid%gdepw_ps(ipt,jpt)) / REAL(irafx)
533          tmp2 = (parentgrid%gdepw_ps(ipt,jpt) - mindepth) / REAL(irafx)
534
535          IF (ABS(slopex) > tmp1) THEN
536             IF (slopex > 0) THEN
537                slopex = tmp1
538             ELSE
539                slopex = -tmp1
540             ENDIF
541          ENDIF
542
543          IF (ABS(slopex) > tmp2) THEN
544             IF (slopex > 0) THEN
545                slopex = tmp2
546             ELSE
547                slopex = -tmp2
548             ENDIF
549          ENDIF
550          !                 
551          ! interpolation on fine grid points (connection zone)
552          !
553          DO ii = i+ii1,i+ii2
554!!             x = ii-i - xdiff/2.
555!!             val = parentgrid%gdepw_ps(ipt,jpt)+slopex * x
556!! chanut: uncomment this to get nearest neighbor interpolation
557             val = parentgrid%gdepw_ps(ipt,jpt)         
558             gdepwtemp(ii,j) = val
559             IF (gdepwtemp(ii,j) < mindepth) THEN
560                gdepwtemp(ii,j) = mindepth
561             ENDIF
562             IF (gdepwtemp(ii,j) > maxdepth) THEN
563                gdepwtemp(ii,j) = maxdepth
564             ENDIF
565             parentbathytab(ii,j) = parentbathy
566          ENDDO
567          ipt =ipt + 1
568       ENDDO
569
570       jpt = jpt + 1               
571    ENDDO
572
573    DO j=jbegin+irafy,jend-irafy,irafy
574
575       DO i=ibegin,iend
576
577          parentbathy = parentbathytab(i,j)
578          IF (parentbathy == 0) THEN
579             mindepth = 0.
580             maxdepth = 0.
581          ELSE
582             mindepth = MAX(gdepw(parentbathy) + MIN( e3zps_min, e3t(parentbathy)*e3zps_rat ),zmin)
583             !                  maxdepth = min(gdepw(parentbathy + 1),zmax)
584             IF (parentbathy < (N-1)) THEN
585                maxdepth = gdepw(parentbathy + 1)
586             ELSE
587                maxdepth = HUGE(1.)
588             ENDIF
589          ENDIF
590
591          slopey = vanleer(gdepwtemp(i,j-irafy:j+irafy:irafy))/REAL(irafy)
592         
593          tmp1 = (maxdepth - gdepwtemp(i,j)) / REAL(irafy)
594          tmp2 = (gdepwtemp(i,j) - mindepth) / REAL(irafy)
595
596          IF (ABS(slopey) > tmp1) THEN
597             IF (slopey > 0) THEN
598                slopey = tmp1
599             ELSE
600                slopey = -tmp1
601             ENDIF
602          ENDIF
603          IF (ABS(slopey) > tmp2) THEN
604             IF (slopey > 0) THEN
605                slopey = tmp2
606             ELSE
607                slopey = -tmp2
608             ENDIF
609          ENDIF
610
611
612          DO jj = j+jj1,j+jj2
613!!             y = jj-j - xdiff/2.
614!!             val = gdepwtemp(i,j) + slopey*y
615!! chanut: uncomment this to get nearest neighbor interpolation
616             val = gdepwtemp(i,j)
617             gdepwtemp(i,jj) = val     
618          ENDDO
619       ENDDO
620    ENDDO
621
622
623    gdepwchild(minboundx:maxboundx,minboundy:maxboundy) = gdepwtemp(minboundx:maxboundx,minboundy:maxboundy)
624    DEALLOCATE(gdepwtemp,parentbathytab)
625
626  END SUBROUTINE correct_level
627  !
628  !
629  !***************************************************
630  ! function van leer to compute the corresponding
631  ! Van Leer slopes
632  !***************************************************
633  !     
634  REAL FUNCTION vanleer(tab)
635    REAL, DIMENSION(3) :: tab
636    REAL res,res1
637    REAL p1,p2,p3
638
639    p1=(tab(3)-tab(1))/2.
640    p2=(tab(2)-tab(1))
641    p3=(tab(3)-tab(2))
642
643    IF ((p1>0.).AND.(p2>0.).AND.(p3>0)) THEN
644       res1=MINVAL((/p1,p2,p3/))
645    ELSEIF ((p1<0.).AND.(p2<0.).AND.(p3<0)) THEN
646       res1=MAXVAL((/p1,p2,p3/))
647    ELSE
648       res1=0.
649    ENDIF
650
651    vanleer = res1   
652
653
654  END FUNCTION vanleer
655  !
656  !
657  !********************************************************************************
658  !   subroutine bathymetry_control                               *
659  !                                                            *
660  !  - Purpose :   check the bathymetry in levels                       *
661  !                              *
662  !  - Method  :   The array mbathy is checked to verified its consistency *
663  !      with the model options. in particular:             *
664  !            mbathy must have at least 1 land grid-points (mbathy<=0)    *
665  !                  along closed boundary.              *
666  !            mbathy must be cyclic IF jperio=1.              *
667  !            mbathy must be lower or equal to jpk-1.            *
668  !            isolated ocean grid points are suppressed from mbathy    *
669  !                  since they are only connected to remaining         *
670  !                  ocean through vertical diffusion.            *
671  !                              *
672  !                              *
673  !********************************************************************************
674
675  SUBROUTINE bathymetry_control(mbathy)
676
677    INTEGER ::   i, j, jl           
678    INTEGER ::   icompt, ibtest, ikmax         
679    REAL*8, DIMENSION(:,:) :: mbathy     
680
681    ! ================
682    ! Bathymetry check
683    ! ================
684
685    ! Suppress isolated ocean grid points
686
687    WRITE(*,*)'                   suppress isolated ocean grid points'
688    WRITE(*,*)'                   -----------------------------------'
689
690    icompt = 0
691
692    DO jl = 1, 2
693       !
694       DO j = 2, SIZE(mbathy,2)-1
695          DO i = 2, SIZE(mbathy,1)-1
696
697             ibtest = MAX( mbathy(i-1,j), mbathy(i+1,j),mbathy(i,j-1),mbathy(i,j+1) )
698             !               
699             IF( ibtest < mbathy(i,j) ) THEN
700                !                 
701                WRITE(*,*) 'grid-point(i,j)= ',i,j,'is changed from',mbathy(i,j),' to ', ibtest
702                mbathy(i,j) = ibtest
703                icompt = icompt + 1
704                !
705             ENDIF
706             !           
707          END DO
708       END DO
709       !
710    END DO
711    !     
712    IF( icompt == 0 ) THEN
713       WRITE(*,*)'     no isolated ocean grid points'
714    ELSE
715       WRITE(*,*)'    ',icompt,' ocean grid points suppressed'
716    ENDIF
717    !
718
719    ! Number of ocean level inferior or equal to jpkm1
720
721    ikmax = 0
722    DO j = 1, SIZE(mbathy,2)
723       DO ji = 1, SIZE(mbathy,1)
724          ikmax = MAX( ikmax, NINT(mbathy(i,j)) )
725       END DO
726    END DO
727    !
728    IF( ikmax > N-1 ) THEN
729       WRITE(*,*) ' maximum number of ocean level = ', ikmax,' >  jpk-1'
730       WRITE(*,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
731    ELSE IF( ikmax < N-1 ) THEN
732       WRITE(*,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 
733       WRITE(*,*) ' you can decrease jpk to ', ikmax+1
734    ENDIF
735
736  END SUBROUTINE bathymetry_control
737  !
738  !
739  !**********************************************************************************
740  !
741  !subroutine get_scale_factors
742  !
743  !**********************************************************************************
744  !
745  SUBROUTINE get_scale_factors(Grid,fse3t,fse3u,fse3v)
746    !
747    IMPLICIT NONE
748    !       
749    TYPE(Coordinates) :: Grid 
750    REAL*8, DIMENSION(:,:,:) :: fse3u,fse3t,fse3v
751    !                                 
752    REAL*8 :: za2,za1,za0,zsur,zacr,zkth,zacr2,zkth2,zdepth,zdepwp,zmin,zmax,zdiff,ze3tp,ze3wp
753    INTEGER :: i,j,jk,jj,ji,jpj,jpi,ik,ii,ipt,jpt,jpk
754    INTEGER, DIMENSION(1) :: k
755    INTEGER :: k1
756    REAL*8, POINTER, DIMENSION(:) :: gdepw,gdept,e3w,e3t
757    REAL*8, POINTER, DIMENSION(:,:)   :: hdepw,e3tp,e3wp
758    REAL*8, POINTER, DIMENSION(:,:,:) :: gdept_ps,gdepw_ps   
759    !       
760    jpi = SIZE(fse3t,1)
761    jpj = SIZE(fse3t,2) 
762    jpk = SIZE(fse3t,3)     
763    !       
764    ALLOCATE(gdepw(jpk),e3t(jpk))
765    ALLOCATE(gdepw_ps(jpi,jpj,jpk))                 
766    ! 
767    IF ( ( pa0 == 0 .OR. pa1 == 0 .OR. psur == 0 ) &
768         .AND. ppdzmin.NE.0 .AND. pphmax.NE.0 ) THEN 
769       !   
770       WRITE(*,*) 'psur,pa0,pa1 computed'
771       za1=( ppdzmin - pphmax / (jpk-1) )          &
772            / ( TANH((1-ppkth)/ppacr) - ppacr/(jpk-1) &
773            *  (  LOG( COSH( (jpk - ppkth) / ppacr) )      &
774            - LOG( COSH( ( 1  - ppkth) / ppacr) )  )  )
775
776       za0  = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
777       zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr )  )
778       !
779    ELSE IF ( (ppdzmin == 0 .OR. pphmax == 0) .AND. psur.NE.0 .AND. &
780         pa0.NE.0 .AND. pa1.NE.0 ) THEN
781       !       
782       WRITE(*,*) 'psur,pa0,pa1 given by namelist'
783       zsur = psur
784       za0  = pa0
785       za1  = pa1
786       za2  = pa2
787       !
788    ELSE
789       !       
790       WRITE(*,*) 'ERROR ***** bad vertical grid parameters ...' 
791       WRITE(*,*) ' '
792       WRITE(*,*) 'please check values of variables'
793       WRITE(*,*) 'in namelist vertical_grid section'
794       WRITE(*,*) ' ' 
795       STOP   
796       !       
797    ENDIF
798
799    zacr  = ppacr
800    zkth  = ppkth     
801    zacr2 = ppacr2
802    zkth2 = ppkth2 
803    !
804    IF( ppkth == 0. ) THEN            !  uniform vertical grid
805         za1 = pphmax / FLOAT(jpk-1) 
806         DO i = 1, jpk
807            gdepw(i) = ( i - 1   ) * za1
808            e3t  (i) =  za1
809         END DO
810    ELSE                            ! Madec & Imbard 1996 function
811       IF( .NOT. ldbletanh ) THEN
812          DO i = 1,jpk
813             !
814             gdepw(i) = (zsur+za0*i+za1*zacr*LOG(COSH((i-zkth)/zacr))) 
815             e3t(i)   = (za0 + za1 * TANH(((i+0.5)-zkth)/zacr))
816             !
817          END DO
818       ELSE
819            DO i = 1,jpk
820               ! Double tanh function
821               gdepw(i) = ( zsur + za0*i  + za1 * zacr * LOG ( COSH( (i-zkth ) / zacr  ) )               &
822                  &                       + za2 * zacr2* LOG ( COSH( (i-zkth2) / zacr2 ) )  )
823               e3t  (i) =          za0         + za1        * TANH(       ((i+0.5)-zkth ) / zacr  )      &
824                  &                            + za2        * TANH(       ((i+0.5)-zkth2) / zacr2 )
825            END DO
826       ENDIF
827    ENDIF         
828    !         
829    gdepw(1)=0.
830    IF ( ln_e3_dep ) THEN      ! e3. = dk[gdep]   
831       !                     
832       DO i = 1, jpk-1
833          e3t(i) = gdepw(i+1)-gdepw(i)
834       END DO
835       e3t(jpk) = e3t(jpk-1)
836    END IF
837    !               
838    DO i = 1,jpk
839       !       
840       fse3t(:,:,i) = e3t(i)
841       gdepw_ps(:,:,i) = gdepw(i)
842       !
843    END DO
844    !                                 
845    gdepw(1) = 0.0 
846    gdepw_ps(:,:,1) = 0.0 
847    !
848    zmax = gdepw(jpk) + e3t(jpk)
849    IF( rn_hmin < 0. ) THEN  ;   i = - INT( rn_hmin )                                  ! from a nb of level
850    ELSE                     ;   i = MINLOC( gdepw, mask = gdepw > rn_hmin, dim = 1 )  ! from a depth
851    ENDIF
852    zmin = gdepw(i+1)
853    !
854    DO jj = 1, jpj
855       DO ji= 1, jpi
856          IF( Grid%bathy_meter(ji,jj) <= 0. ) THEN
857             Grid%bathy_meter(ji,jj) = 0.e0
858          ELSE
859             Grid%bathy_meter(ji,jj) = MAX( Grid%bathy_meter(ji,jj), zmin )
860             Grid%bathy_meter(ji,jj) = MIN( Grid%bathy_meter(ji,jj), zmax )
861          ENDIF
862       END DO
863    END DO
864    !
865    DO jj = 1, jpj
866       DO ji = 1, jpi
867          ik = Grid%bathy_level(ji,jj) 
868          IF( ik > 0 ) THEN
869             ! max ocean level case
870             IF( ik == jpk-1 ) THEN
871                zdepwp = Grid%bathy_meter(ji,jj)
872                ze3tp  = Grid%bathy_meter(ji,jj) - gdepw(ik)
873                fse3t(ji,jj,ik  ) = ze3tp
874                fse3t(ji,jj,ik+1) = ze3tp
875                gdepw_ps(ji,jj,ik+1) = zdepwp
876             ELSE
877                IF( Grid%bathy_meter(ji,jj) <= gdepw(ik+1) ) THEN
878                   gdepw_ps(ji,jj,ik+1) = Grid%bathy_meter(ji,jj)
879                ELSE
880                   gdepw_ps(ji,jj,ik+1) = gdepw(ik+1)
881                ENDIF
882                fse3t(ji,jj,ik) = e3t(ik) * ( gdepw_ps(ji,jj,ik+1) - gdepw(ik))   & 
883                     /( gdepw(ik+1) - gdepw(ik)) 
884                fse3t(ji,jj,ik+1) = fse3t(ji,jj,ik)
885
886             ENDIF
887          ENDIF
888       END DO
889    END DO
890    !
891    DO i = 1, jpk
892       fse3u (:,:,i)  = e3t(i)
893       fse3v (:,:,i)  = e3t(i)
894    END DO
895    !
896    DO jk = 1,jpk
897       DO jj = 1, jpj-1
898          DO ji = 1, jpi-1
899             fse3u (ji,jj,jk) = MIN( fse3t(ji,jj,jk), fse3t(ji+1,jj,jk))
900             fse3v (ji,jj,jk) = MIN( fse3t(ji,jj,jk), fse3t(ji,jj+1,jk))     
901          ENDDO
902       ENDDO
903    ENDDO
904    !           
905    DEALLOCATE(gdepw,e3t)
906    DEALLOCATE(gdepw_ps)
907    DEALLOCATE(Grid%bathy_meter,Grid%bathy_level) 
908    !
909  END SUBROUTINE get_scale_factors
910  !       
911END MODULE agrif_partial_steps
912
913
Note: See TracBrowser for help on using the repository browser.