Changeset 13204 for utils/tools/DOMAINcfg/src/domzgr.F90
- Timestamp:
- 2020-07-02T10:38:35+02:00 (4 years ago)
- File:
-
- 1 edited
-
utils/tools/DOMAINcfg/src/domzgr.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
utils/tools/DOMAINcfg/src/domzgr.F90
r12414 r13204 17 17 !! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function 18 18 !! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capability 19 !! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye 20 20 !! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying 21 21 !!---------------------------------------------------------------------- … … 36 36 !!--------------------------------------------------------------------- 37 37 USE dom_oce ! ocean domain 38 USE depth_e3 ! depth <=> e3 38 39 ! 39 40 USE in_out_manager ! I/O manager … … 44 45 USE dombat 45 46 USE domisf 47 USE agrif_domzgr 46 48 47 49 IMPLICIT NONE … … 92 94 CONTAINS 93 95 94 SUBROUTINE dom_zgr 96 SUBROUTINE dom_zgr( k_top, k_bot ) 95 97 !!---------------------------------------------------------------------- 96 98 !! *** ROUTINE dom_zgr *** … … 109 111 !! ** Action : define gdep., e3., mbathy and bathy 110 112 !!---------------------------------------------------------------------- 113 INTEGER, DIMENSION(:,:), INTENT(out) :: k_top, k_bot ! ocean first and last level indices 114 ! 111 115 INTEGER :: ioptio, ibat ! local integer 112 116 INTEGER :: ios 113 117 ! 114 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 118 INTEGER :: jk 119 REAL(wp) :: zrefdep ! depth of the reference level (~10m) 120 121 122 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 115 123 !!---------------------------------------------------------------------- 116 124 ! … … 124 132 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp ) 125 133 IF(lwm) WRITE ( numond, namzgr ) 134 135 IF(ln_read_cfg) THEN 136 IF(lwp) WRITE(numout,*) 137 IF(lwp) WRITE(numout,*) ' ==>>> Read vertical mesh in ', TRIM( cn_domcfg ), ' file' 138 ! 139 CALL zgr_read ( ln_zco , ln_zps , ln_sco, ln_isfcav, & 140 & gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth 141 & gdept_0 , gdepw_0 , & ! gridpoints depth 142 & e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors 143 & e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors 144 & k_top , k_bot ) ! 1st & last ocean level 145 ! 146 !!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears 147 ! ! Compute gde3w_0 (vertical sum of e3w) 148 ! gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 149 ! DO jk = 2, jpk 150 ! gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 151 ! END DO 152 ! 153 154 ! ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top) 155 CALL zgr_top_bot( k_top, k_bot ) ! with a minimum value set to 1 156 157 ! ! deepest/shallowest W level Above/Below ~10m 158 !!gm BUG in s-coordinate this does not work! 159 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness) 160 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 161 nla10 = nlb10 - 1 ! deepest W level Above ~10m 162 !!gm end bug 163 164 ENDIF 126 165 127 166 IF(lwp) THEN ! Control print … … 134 173 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 135 174 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 136 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 137 ENDIF 138 139 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 175 ENDIF 140 176 141 177 ioptio = 0 ! Check Vertical coordinate options … … 150 186 IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1 151 187 IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' ) 188 189 IF(.NOT.ln_read_cfg) THEN 152 190 ! 153 191 ! Build the vertical coordinate system … … 163 201 ! ----------------------------------- 164 202 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 203 k_bot = mbkt 165 204 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 166 ! 167 IF( nprint == 1 .AND. lwp ) THEN 168 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 205 k_top = mikt 206 WHERE( bathy(:,:) <= 0._wp ) k_top(:,:) = 0 ! set k_top to zero over land 207 ENDIF 208 ! 209 IF( lwp ) THEN 210 WRITE(numout,*) ' MIN val k_top ', MINVAL( k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) ) 211 WRITE(numout,*) ' MIN val k_bot ', MINVAL( k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) ) 169 212 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 170 213 & ' w ', MINVAL( gdepw_0(:,:,:) ) … … 181 224 & ' w ', MAXVAL( e3w_0(:,:,:) ) 182 225 ENDIF 226 183 227 ! 184 228 END SUBROUTINE dom_zgr 185 229 230 SUBROUTINE zgr_read( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate 231 & pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate 232 & pdept , pdepw , & ! 3D t & w-points depth 233 & pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors 234 & pe3w , pe3uw , pe3vw , & ! - - - 235 & k_top , k_bot ) ! top & bottom ocean level 236 !!--------------------------------------------------------------------- 237 !! *** ROUTINE zgr_read *** 238 !! 239 !! ** Purpose : Read the vertical information in the domain configuration file 240 !! 241 !!---------------------------------------------------------------------- 242 LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags 243 LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag 244 REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m] 245 REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m] 246 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m] 247 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m] 248 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - - 249 INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level 250 ! 251 INTEGER :: jk ! dummy loop index 252 INTEGER :: inum ! local logical unit 253 REAL(WP) :: z_zco, z_zps, z_sco, z_cav 254 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 255 !!---------------------------------------------------------------------- 256 ! 257 IF(lwp) THEN 258 WRITE(numout,*) 259 WRITE(numout,*) ' zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file' 260 WRITE(numout,*) ' ~~~~~~~~' 261 ENDIF 262 ! 263 CALL iom_open( cn_domcfg, inum ) 264 ! 265 ! !* type of vertical coordinate 266 CALL iom_get( inum, 'ln_zco' , z_zco ) 267 CALL iom_get( inum, 'ln_zps' , z_zps ) 268 CALL iom_get( inum, 'ln_sco' , z_sco ) 269 IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF 270 IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF 271 IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF 272 ! 273 ! !* ocean cavities under iceshelves 274 CALL iom_get( inum, 'ln_isfcav', z_cav ) 275 IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF 276 ! 277 ! !* vertical scale factors 278 CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate 279 CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d ) 280 ! 281 CALL iom_get( inum, jpdom_data, 'e3t_0' , pe3t , lrowattr=ln_use_jattr ) ! 3D coordinate 282 CALL iom_get( inum, jpdom_data, 'e3u_0' , pe3u , lrowattr=ln_use_jattr ) 283 CALL iom_get( inum, jpdom_data, 'e3v_0' , pe3v , lrowattr=ln_use_jattr ) 284 CALL iom_get( inum, jpdom_data, 'e3f_0' , pe3f , lrowattr=ln_use_jattr ) 285 CALL iom_get( inum, jpdom_data, 'e3w_0' , pe3w , lrowattr=ln_use_jattr ) 286 CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw , lrowattr=ln_use_jattr ) 287 CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw , lrowattr=ln_use_jattr ) 288 ! 289 ! !* depths 290 ! !- old depth definition (obsolescent feature) 291 IF( iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0 .AND. & 292 & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0 .AND. & 293 & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0 .AND. & 294 & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0 ) THEN 295 CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 296 & ' depths at t- and w-points read in the domain configuration file') 297 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d ) 298 CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d ) 299 CALL iom_get( inum, jpdom_data , 'gdept_0' , pdept , lrowattr=ln_use_jattr ) 300 CALL iom_get( inum, jpdom_data , 'gdepw_0' , pdepw , lrowattr=ln_use_jattr ) 301 ! 302 ELSE !- depths computed from e3. scale factors 303 CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) ! 1D reference depth 304 CALL e3_to_depth( pe3t , pe3w , pdept , pdepw ) ! 3D depths 305 IF(lwp) THEN 306 WRITE(numout,*) 307 WRITE(numout,*) ' Reference 1D z-coordinate depth and scale factors:' 308 WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" ) 309 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk ) 310 ENDIF 311 ENDIF 312 ! 313 ! !* ocean top and bottom level 314 CALL iom_get( inum, jpdom_data, 'top_level' , z2d , lrowattr=ln_use_jattr ) ! 1st wet T-points (ISF) 315 k_top(:,:) = NINT( z2d(:,:) ) 316 CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d , lrowattr=ln_use_jattr ) ! last wet T-points 317 k_bot(:,:) = NINT( z2d(:,:) ) 318 ! 319 ! reference depth for negative bathy (wetting and drying only) 320 ! IF( ll_wd ) CALL iom_get( inum, 'rn_wd_ref_depth' , ssh_ref ) 321 ! 322 CALL iom_close( inum ) 323 ! 324 END SUBROUTINE zgr_read 325 326 SUBROUTINE zgr_top_bot( k_top, k_bot ) 327 !!---------------------------------------------------------------------- 328 !! *** ROUTINE zgr_top_bot *** 329 !! 330 !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) 331 !! 332 !! ** Method : computes from k_top and k_bot with a minimum value of 1 over land 333 !! 334 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 335 !! ocean level at t-, u- & v-points 336 !! (min value = 1) 337 !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest 338 !! ocean level at t-, u- & v-points 339 !! (min value = 1 over land) 340 !!---------------------------------------------------------------------- 341 INTEGER , DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! top & bottom ocean level indices 342 ! 343 INTEGER :: ji, jj ! dummy loop indices 344 REAL(wp), DIMENSION(jpi,jpj) :: zk ! workspace 345 !!---------------------------------------------------------------------- 346 ! 347 IF(lwp) WRITE(numout,*) 348 IF(lwp) WRITE(numout,*) ' zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels ' 349 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 350 ! 351 mikt(:,:) = MAX( k_top(:,:) , 1 ) ! top ocean k-index of T-level (=1 over land) 352 ! 353 mbkt(:,:) = MAX( k_bot(:,:) , 1 ) ! bottom ocean k-index of T-level (=1 over land) 354 355 ! ! N.B. top k-index of W-level = mikt 356 ! ! bottom k-index of W-level = mbkt+1 357 DO jj = 1, jpjm1 358 DO ji = 1, jpim1 359 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 360 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 361 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 362 ! 363 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) 364 mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) 365 END DO 366 END DO 367 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 368 zk(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; miku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 369 zk(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mikv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 370 zk(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'F', 1. ) ; mikf(:,:) = MAX( NINT( zk(:,:) ), 1 ) 371 ! 372 zk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; mbku(:,:) = MAX( NINT( zk(:,:) ), 1 ) 373 zk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 ) 374 ! 375 END SUBROUTINE zgr_top_bot 186 376 187 377 SUBROUTINE zgr_z … … 665 855 ENDIF 666 856 667 zbathy(:,:) = FLOAT( mbathy(:,:) ) 668 CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp ) 669 mbathy(:,:) = INT( zbathy(:,:) ) 670 857 IF( lk_mpp ) THEN 858 zbathy(:,:) = FLOAT( mbathy(:,:) ) 859 CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp ) 860 mbathy(:,:) = INT( zbathy(:,:) ) 861 ENDIF 671 862 ! ! East-west cyclic boundary conditions 672 863 IF( jperio == 0 ) THEN … … 1061 1252 END SUBROUTINE zgr_zps 1062 1253 1254 1063 1255 SUBROUTINE zgr_sco 1064 1256 !!---------------------------------------------------------------------- … … 1262 1454 END DO 1263 1455 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1264 CALL lbc_lnk( ' domzgr',zenv, 'T', 1._wp, 'no0' )1456 CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' ) 1265 1457 ! ! ================ ! 1266 1458 END DO ! End loop ! … … 1341 1533 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 1342 1534 1343 IF( nprint == 1 .AND.lwp ) THEN1535 IF( lwp ) THEN 1344 1536 WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & 1345 1537 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) … … 1397 1589 END DO 1398 1590 END DO 1399 IF( nprint == 1 .AND.lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &1591 IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 1400 1592 & ' MAX ', MAXVAL( mbathy(:,:) ) 1401 1593 1402 IF( nprint == 1 .AND.lwp ) THEN ! min max values over the local domain1594 IF( lwp ) THEN ! min max values over the local domain 1403 1595 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1404 1596 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & … … 1788 1980 z_gsigt(jk) = -fssig( REAL(jk,wp) ) 1789 1981 END DO 1790 IF( nprint == 1 .AND.lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)1982 IF( lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk) 1791 1983 ! 1792 1984 ! Coefficients for vertical scale factors at w-, t- levels
Note: See TracChangeset
for help on using the changeset viewer.
