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.
Changeset 15279 for utils/tools/DOMAINcfg/src/agrif_connect.F90 – NEMO

Ignore:
Timestamp:
2021-09-23T12:00:23+02:00 (3 years ago)
Author:
jchanut
Message:

#2222 and #2638: Enable creating agrif meshes with different vertical grids (geopotential only as a start)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • utils/tools/DOMAINcfg/src/agrif_connect.F90

    r14952 r15279  
    22 
    33   USE dom_oce 
    4    USE domzgr 
    54   USE agrif_parameters 
    65   USE agrif_profiles 
     
    98   PRIVATE 
    109 
    11    PUBLIC agrif_boundary_connections  
     10   PUBLIC agrif_boundary_connections, agrif_bathymetry_connect  
    1211 
    1312CONTAINS 
     
    2726!      CALL Agrif_Bc_variable(e3t_copy_id, procname = connect_e3t_copy) 
    2827 
    29       ALLOCATE(e3t_interp(jpi,jpj,jpk)) 
    30       e3t_interp = -10. 
     28      ALLOCATE(e3t_interp_done(jpi,jpj)) 
     29      e3t_interp_done(:,:) = .FALSE.  
     30      ! set extrapolation on for interpolation near the coastline: 
     31      Agrif_UseSpecialValue = .TRUE. 
     32      Agrif_SpecialValue = 0._wp 
     33      CALL Agrif_Bc_variable(e3t_connect_id, procname = connect_e3t_connect) 
     34      ! Override in ghost zone by nearest value: 
    3135      Agrif_UseSpecialValue = .FALSE. 
    32       Agrif_SpecialValue = 0. 
    33       CALL Agrif_Bc_variable(e3t_connect_id, procname = connect_e3t_connect) 
     36      e3t_interp_done(:,:) = .FALSE. 
     37      CALL Agrif_Bc_variable(e3t_copy_id,    procname = connect_e3t_connect) 
    3438      Agrif_UseSpecialValue = .FALSE. 
     39      DEALLOCATE(e3t_interp_done) 
    3540      ! 
    3641   END SUBROUTINE agrif_boundary_connections 
     42 
     43   SUBROUTINE agrif_bathymetry_connect 
     44      !!---------------------------------------------------------------------- 
     45      !!                  ***  ROUTINE agrif_bathymetry_connect  *** 
     46      !!----------------------------------------------------------------------   
     47      IF( Agrif_Root() ) return 
     48 
     49      CALL agrif_connection() 
     50      ! 
     51      ALLOCATE(e3t_interp_done(jpi,jpj)) 
     52      e3t_interp_done(:,:) = .FALSE.  
     53      ! set extrapolation on for interpolation near the coastline: 
     54      Agrif_UseSpecialValue = .TRUE. 
     55      Agrif_SpecialValue = 0._wp 
     56      CALL Agrif_Bc_variable(e3t_connect_id, procname = connect_bathy_connect) 
     57      ! Override in ghost zone by nearest value: 
     58      Agrif_UseSpecialValue = .FALSE. 
     59      e3t_interp_done(:,:) = .FALSE. 
     60      CALL Agrif_Bc_variable(e3t_copy_id,    procname = connect_bathy_connect) 
     61      Agrif_UseSpecialValue = .FALSE. 
     62      DEALLOCATE(e3t_interp_done) 
     63      ! 
     64   END SUBROUTINE agrif_bathymetry_connect 
    3765 
    3866   SUBROUTINE connect_e3t_copy( ptab, i1, i2, j1, j2, k1, k2, before, nb,ndir) 
     
    5078         ptab(i1:i2,j1:j2,k1:k2) = e3t_0(i1:i2,j1:j2,k1:k2) 
    5179      ELSE 
    52          e3t_0(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) 
     80         e3t_0(i1:i2,j1:j2,1:jpk) = ptab(i1:i2,j1:j2,1:jpk) 
    5381      ENDIF 
    5482      ! 
     
    89117      ! 
    90118      !!----------------------------------------------------------------------  
    91       INTEGER :: ji, jj, jk  
     119      INTEGER :: ji, jj, jk, ik  
    92120      REAL(wp), DIMENSION(i1:i2,j1:j2) :: bathy_local, bathy_interp 
    93       REAL(wp) :: zdepth, zmax  
     121      REAL(wp) :: zdepth, zdepwp, zmax, ze3tp, ze3wp, zhmin  
    94122      ! 
    95123      IF( before) THEN 
    96          DO jk=1,jpk 
     124         DO jk=k1, k2 
    97125            DO jj=j1,j2 
    98126               DO ji=i1,i2 
     
    108136         DO jj=j1,j2 
    109137            DO ji=i1,i2 
    110                ptab(ji,jj,jpk+1) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 
     138               ptab(ji,jj,k2) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 
    111139            END DO 
    112140         END DO 
     
    115143            DO ji=i1,i2 
    116144               bathy_local (ji,jj) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 
    117                bathy_interp (ji,jj) = ptab(ji,jj,jpk+1) 
    118  
     145               bathy_interp (ji,jj) = ptab(ji,jj,k2) 
     146               ! keep child masking in transition zone: 
     147               IF ((ztabramp(ji,jj)/=1._wp).AND.(bathy_local(ji,jj)==0._wp)) bathy_interp(ji,jj)=0._wp 
    119148        ! Connected bathymetry 
    120                IF( e3t_interp(ji,jj,1) == -10 ) THEN 
     149               IF( .NOT.e3t_interp_done(ji,jj) ) THEN 
    121150                  bathy_local(ji,jj)=(1.-ztabramp(ji,jj))*bathy_local(ji,jj)+ztabramp(ji,jj)*bathy_interp(ji,jj) 
    122151               ENDIF 
     
    125154 
    126155        ! Update mbkt and ssmask 
     156         IF( rn_hmin < 0._wp ) THEN 
     157            ik = - INT( rn_hmin ) 
     158         ELSE 
     159            ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) 
     160         ENDIF 
     161         zhmin = gdepw_1d(ik+1) 
     162 
    127163         zmax = gdepw_1d(jpk) + e3t_1d(jpk) 
    128164         bathy_local(:,:) = MAX(MIN(zmax,bathy_local(:,:)),0._wp) 
    129          WHERE( bathy_local(i1:i2,j1:j2) == 0._wp); mbathy(i1:i2,j1:j2) = 0 
    130          ELSE WHERE                       ; mbathy(i1:i2,j1:j2) = jpkm1 
     165         WHERE( bathy_local(i1:i2,j1:j2) == 0._wp) 
     166            mbathy(i1:i2,j1:j2) = 0 
     167         ELSE WHERE  
     168            mbathy(i1:i2,j1:j2) = jpkm1  
     169            bathy_local(i1:i2,j1:j2) = MAX(  zhmin , bathy_local(i1:i2,j1:j2)  )  
    131170         END WHERE 
    132171 
    133172         DO jk=jpkm1,1,-1 
    134            zdepth = gdepw_1d(jk)+MIN(e3zps_min,e3t_1d(jk)*e3zps_rat) 
     173           zdepth = gdepw_1d(jk) + MIN(e3zps_min,e3t_1d(jk)*e3zps_rat) 
    135174           WHERE( 0._wp < bathy_local(:,:) .AND. bathy_local(:,:) <= zdepth ) mbathy(i1:i2,j1:j2) = jk-1 
    136175         ENDDO 
     
    141180          
    142181         mbkt(i1:i2,j1:j2) = MAX( mbathy(i1:i2,j1:j2), 1 ) 
    143  
    144182         ! 
    145          DO jk=1,jpk 
     183         DO jj = j1, j2 
     184            DO ji = i1, i2 
     185               IF( .NOT.e3t_interp_done(ji,jj) ) THEN ! the connection has not yet been done 
     186                  DO jk = 1, jpk     
     187                     gdept_0(ji,jj,jk) = gdept_1d(jk) 
     188                     gdepw_0(ji,jj,jk) = gdepw_1d(jk) 
     189                     e3t_0  (ji,jj,jk) = e3t_1d  (jk) 
     190                     e3w_0  (ji,jj,jk) = e3w_1d  (jk) 
     191                  END DO  
     192                  ! 
     193                  ik = mbathy(ji,jj) 
     194                  IF( ik > 0 ) THEN               ! ocean point only 
     195                     ! max ocean level case 
     196                     IF( ik == jpkm1 ) THEN 
     197                        zdepwp = bathy_local(ji,jj) 
     198                        ze3tp  = bathy_local(ji,jj) - gdepw_1d(ik) 
     199                        ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 
     200                        e3t_0(ji,jj,ik  ) = ze3tp 
     201                        e3t_0(ji,jj,ik+1) = ze3tp 
     202                        e3w_0(ji,jj,ik  ) = ze3wp 
     203                        e3w_0(ji,jj,ik+1) = ze3tp 
     204                        gdepw_0(ji,jj,ik+1) = zdepwp 
     205                        gdept_0(ji,jj,ik  ) = gdept_1d(ik-1) + ze3wp 
     206                        gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 
     207                        ! 
     208                     ELSE                         ! standard case 
     209                        IF( bathy_local(ji,jj) <= gdepw_1d(ik+1) ) THEN 
     210                           gdepw_0(ji,jj,ik+1) = bathy_local(ji,jj) 
     211                        ELSE 
     212                           gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 
     213                        ENDIF 
     214                        gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 
     215                              &                * ((gdept_1d(     ik  ) - gdepw_1d(ik) )           & 
     216                              &                / ( gdepw_1d(     ik+1) - gdepw_1d(ik) )) 
     217                        e3t_0  (ji,jj,ik) = e3t_1d  (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik)) & 
     218                              &                / ( gdepw_1d(      ik+1) - gdepw_1d(ik)) 
     219                        e3w_0(ji,jj,ik) = & 
     220                              & 0.5_wp * (gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) )   & 
     221                              &        * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 
     222                        !       ... on ik+1 
     223                        e3w_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     224                        e3t_0  (ji,jj,ik+1) = e3t_0  (ji,jj,ik) 
     225                        gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 
     226                     ENDIF 
     227                  ENDIF 
     228               ENDIF 
     229               e3t_interp_done(ji,jj) = .TRUE. 
     230            END DO 
     231         END DO 
     232      ENDIF 
     233      ! 
     234   END SUBROUTINE connect_e3t_connect 
     235 
     236   SUBROUTINE connect_bathy_connect( ptab, i1, i2, j1, j2, k1, k2, before, nb,ndir) 
     237      !!---------------------------------------------------------------------- 
     238      !!                  ***  ROUTINE connect_e3t_connect  *** 
     239      !!----------------------------------------------------------------------   
     240      INTEGER                               , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2 
     241      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) ::   ptab 
     242      LOGICAL                               , INTENT(in   ) ::   before 
     243      INTEGER                               , INTENT(in   ) ::   nb , ndir 
     244      ! 
     245      !!----------------------------------------------------------------------  
     246      INTEGER :: ji, jj, jk 
     247      ! 
     248      IF( before) THEN 
     249         DO jk=k1,k2 
    146250            DO jj=j1,j2 
    147251               DO ji=i1,i2 
    148                   IF( e3t_interp(ji,jj,jk) == -10 ) THEN ! the connection has not yet been done 
    149                      e3t_interp(ji,jj,jk) = MAX( ptab(ji,jj,jk),MIN(e3zps_min, e3t_1d(jk)*e3zps_rat) ) 
    150                   !   e3t_interp(ji,jj,jk) = MIN( e3t_interp(ji,jj,jk),e3t_1d(jk) ) 
    151                      e3t_0(ji,jj,jk) = ( 1. - ztabramp(ji,jj) )*e3t_0(ji,jj,jk) + ztabramp(ji,jj)*e3t_interp(ji,jj,jk) 
     252                  IF( mbkt(ji,jj) .GE. jk ) THEN 
     253                     ptab(ji,jj,jk) = e3t_0(ji,jj,jk) 
     254                  ELSE 
     255                     ptab(ji,jj,jk) = 0._wp 
    152256                  ENDIF 
    153                   IF( jk > mbkt(ji,jj)) THEN 
    154                     e3t_0(ji,jj,jk) = e3t_1d(jk) 
    155                   ENDIF 
    156              END DO 
    157            END DO 
    158          END DO 
    159       ENDIF 
    160       ! 
    161    END SUBROUTINE connect_e3t_connect 
     257               END DO 
     258            END DO 
     259         END DO 
     260         ! 
     261         DO jj=j1,j2 
     262            DO ji=i1,i2 
     263               ptab(ji,jj,k2) = SUM ( e3t_0(ji,jj, 1:mbkt(ji,jj) ) ) * ssmask(ji,jj) 
     264            END DO 
     265         END DO 
     266      ELSE 
     267         DO jj=j1,j2 
     268            DO ji=i1,i2 
     269               ! keep child masking in transition zone: 
     270               IF ((ztabramp(ji,jj)/=1._wp).AND.(bathy(ji,jj)==0._wp)) ptab(ji,jj,k2) = 0._wp 
     271               ! Connected bathymetry 
     272               IF( .NOT.e3t_interp_done(ji,jj) ) THEN 
     273                  bathy(ji,jj)=(1._wp-ztabramp(ji,jj))*bathy(ji,jj)+ztabramp(ji,jj)*ptab(ji,jj,k2) 
     274                  e3t_interp_done(ji,jj) = .TRUE. 
     275               ENDIF 
     276            END DO 
     277         END DO 
     278      ENDIF 
     279      ! 
     280   END SUBROUTINE connect_bathy_connect 
    162281    
    163282   SUBROUTINE agrif_connection 
     
    181300 
    182301      ! --- West --- ! 
    183       IF( ((nbondi == -1) .OR. (nbondi == 2) ).AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6)) THEN 
     302      IF( lk_west ) THEN 
    184303         ind1 = nn_hls + nbghostcells + istart 
    185304         ind2 = ind1 + ispongearea  
     
    200319 
    201320      ! --- East --- ! 
    202       IF( ((nbondi == 1) .OR. (nbondi == 2) ).AND. .NOT. (jperio == 1 .OR. jperio == 4 .OR. jperio == 6)) THEN 
     321      IF( lk_east ) THEN 
    203322         ind2 = jpiglo -  (nn_hls + nbghostcells -1 ) - istart 
    204323         ind1 = ind2 -ispongearea        
     
    223342 
    224343      ! --- South --- ! 
    225       IF(( (nbondj == -1) .OR. (nbondj == 2) ).AND.(lk_south)) THEN 
     344      IF( lk_south ) THEN 
    226345         ind1 = nn_hls + nbghostcells + istart 
    227346         ind2 = ind1 + ispongearea  
     
    242361 
    243362      ! --- North --- ! 
    244       IF(( (nbondj == 1) .OR. (nbondj == 2) ).AND.(lk_north)) THEN 
     363      IF( lk_north ) THEN 
    245364         ind2 = jpjglo - (nn_hls + nbghostcells - 1) - istart 
    246365         ind1 = ind2 -ispongearea 
     
    265384   SUBROUTINE agrif_boundary_connections 
    266385   END SUBROUTINE agrif_boundary_connections 
     386   SUBROUTINE agrif_bathymetry_connect  
     387   END SUBROUTINE agrif_bathymetry_connect  
    267388#endif 
    268389 
Note: See TracChangeset for help on using the changeset viewer.