Changeset 15279 for utils/tools/DOMAINcfg/src/domzgr.F90

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:

• utils/tools/DOMAINcfg/src/domzgr.F90 (modified) (8 diffs)

utils/tools/DOMAINcfg/src/domzgr.F90

@@ -45 +45 @@
    USE dombat
    USE domisf
+   USE agrif_connect
    USE agrif_domzgr
 
@@ -176 +177 @@
       IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
       IF( ioptio > 0 )   CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
+
+#if defined key_agrif
+!      IF ( (.NOT.Agrif_Root()).AND.((.NOT.ln_zps).OR.(.NOT.Agrif_parent(ln_zps)))) THEN
+!         CALL ctl_stop( 'STOP', 'AGRIF zooms require ln_zps=T for both Child and Parent')
+!      ENDIF
+#endif
 
       IF(.NOT.ln_read_cfg) THEN
@@ -560 +567 @@
       REAL(wp) ::   r_bump , h_bump , h_oce   ! bump characteristics 
       REAL(wp) ::   zi, zj, zh, zhmin         ! local scalars
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zrand 
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zrand, zbatv 
       !!----------------------------------------------------------------------
       !
@@ -609 +616 @@
             ENDIF
             IF ( cp_cfg=='DOME' ) THEN
-               bathy(:,:) = MIN(3600._wp, MAX( 600._wp,  600._wp -gphit(:,:)*1.e3*0.01 ))
+               ALLOCATE(zbatv(jpi,jpj))
+               zbatv(:,:) = MIN(3600._wp, MAX( 600._wp,  600._wp -gphiv(:,:)*1.e3*0.01 ))
+               bathy(:,1) = 0._wp
+               DO jj =2,jpj
+                  bathy(:,jj) = 0.5_wp*(zbatv(:,jj) + zbatv(:,jj-1))
+               END DO 
+               CALL lbc_lnk( 'zgr_bat', bathy, 'T', 1._wp )
                WHERE (gphit(:,:) >0._wp) bathy(:,:) = 0._wp
                ! Dig inlet:
                WHERE ((gphit(:,:)>0._wp).AND.(glamt(:,:)>-50._wp).AND.(glamt(:,:)<50._wp)) bathy(:,:) = 600._wp
+               DEALLOCATE(zbatv)
             ENDIF
 !            CALL lbc_lnk( 'zgr_bat', bathy, 'T', 1._wp )
@@ -731 +745 @@
       ENDIF
       !
+#if defined key_agrif
+      IF ( .NOT.Agrif_Root() ) CALL agrif_bathymetry_connect
+#endif
+      !
       IF ( .not. ln_sco ) THEN                                !==  set a minimum depth  ==!
          IF( rn_hmin < 0._wp ) THEN    ;   ik = - INT( rn_hmin )                                      ! from a nb of level
@@ -827 +845 @@
            mbathy(  mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = jpkm1
          ENDIF
-      ELSEIF( l_Iperio ) THEN
+      ELSEIF( jperio == 1 .OR. jperio == 4 .OR. jperio ==  6 ) THEN
          IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
-         mbathy( 1 ,:) = mbathy(jpim1,:)
-         mbathy(jpi,:) = mbathy(  2  ,:)
+!         mbathy( 1 ,:) = mbathy(jpim1,:)
+!         mbathy(jpi,:) = mbathy(  2  ,:)
       ELSEIF( jperio == 2 ) THEN
          IF(lwp) WRITE(numout,*) '   equatorial boundary conditions on mbathy: jperio = ', jperio
@@ -1409 +1427 @@
       ! Envelope bathymetry saved in hbatt
       hbatt(:,:) = zenv(:,:) 
-      IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
+      IF ((ntopo>0).AND.MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
          CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
          DO jj = 1, jpj
@@ -1532 +1550 @@
          END DO
       END DO
+
+      WHERE (bathy(:,:)<=0) mbathy(:,:) = 0
+
       IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ),   &
          &                                                       ' MAX ', MAXVAL( mbathy(:,:) )