source: branches/DEV_r1879_FCM/NEMOGCM/TOOLS/NESTING/src/agrif_partial_steps.f90 @ 2143

Last change on this file since 2143 was 2143, checked in by rblod, 11 years ago

Improvement of FCM branch

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