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 @ 12264

Last change on this file since 12264 was 12264, checked in by clem, 5 years ago

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