Changeset 1099 for trunk/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2008-06-09T10:58:04+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (53 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90
r1083 r1099 4 4 !! Ocean initialization : domain initialization 5 5 !!============================================================================== 6 6 !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate 7 !! ! 1997-07 (G. Madec) lbc_lnk call 8 !! ! 1997-04 (J.-O. Beismann) 9 !! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module 10 !! - ! 2002-09 (A. de Miranda) rigid-lid + islands 11 !! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module 12 !! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function 13 !! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco 14 !! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style 7 15 !!---------------------------------------------------------------------- 8 !! dom_zgr : defined the ocean vertical coordinate system 16 17 !!---------------------------------------------------------------------- 18 !! dom_zgr : defined the ocean vertical coordinate system 9 19 !! zgr_bat : bathymetry fields (levels and meters) 10 20 !! zgr_bat_zoom : modify the bathymetry field if zoom domain … … 15 25 !! zgr_sco : s-coordinate 16 26 !! fssig : sigma coordinate non-dimensional function 27 !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!) 17 28 !!--------------------------------------------------------------------- 18 !! * Modules used19 29 USE oce ! ocean dynamics and tracers 20 30 USE dom_oce ! ocean space and time domain … … 22 32 USE lib_mpp ! distributed memory computing library 23 33 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE closea 25 USE solisl 34 USE closea ! closed seas 35 USE solisl ! solver - island in rigid-lid 26 36 USE c1d ! 1D configuration 27 37 … … 29 39 PRIVATE 30 40 31 !! * Routine accessibility 32 PUBLIC dom_zgr ! called by dom_init.F90 33 34 !! * Module variables 35 REAL(wp) :: & !!: Namelist nam_zgr_sco 36 sbot_min = 300. , & !: minimum depth of s-bottom surface (>0) (m) 37 sbot_max = 5250. , & !: maximum depth of s-bottom surface (= ocean depth) (>0) (m) 38 theta = 6.0 , & !: surface control parameter (0<=theta<=20) 39 thetb = 0.75, & !: bottom control parameter (0<=thetb<= 1) 40 r_max = 0.15 !: maximum cut-off r-value allowed (0<r_max<1) 41 42 41 PUBLIC dom_zgr ! called by dom_init.F90 42 43 !!gm DOCTOR name for the namelist parameter... 44 ! !!! ** Namelist nam_zgr_sco ** 45 REAL(wp) :: sbot_min = 300. ! minimum depth of s-bottom surface (>0) (m) 46 REAL(wp) :: sbot_max = 5250. ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 47 REAL(wp) :: theta = 6.0 ! surface control parameter (0<=theta<=20) 48 REAL(wp) :: thetb = 0.75 ! bottom control parameter (0<=thetb<= 1) 49 REAL(wp) :: r_max = 0.15 ! maximum cut-off r-value allowed (0<r_max<1) 43 50 44 51 !! * Substitutions … … 46 53 # include "vectopt_loop_substitute.h90" 47 54 !!---------------------------------------------------------------------- 48 !! OPA 9.0 , LOCEAN-IPSL (2005) 55 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 56 !! $Id:$ 57 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 49 58 !!---------------------------------------------------------------------- 50 59 … … 58 67 !! vertical scale factors. 59 68 !! 60 !! ** Method reference vertical coordinate 61 !! 62 !! ** Action : 63 !! 64 !! History : 65 !! 9.0 ! 03-08 (G. Madec) original code 66 !! 9.0 ! 05-10 (A. Beckmann) modifications for hybrid s-ccordinates 67 !!---------------------------------------------------------------------- 68 INTEGER :: ioptio = 0 ! temporary integer 69 69 !! ** Method : - reference 1D vertical coordinate (gdep._0, e3._0) 70 !! - read/set ocean depth and ocean levels (bathy, mbathy) 71 !! - vertical coordinate (gdep., e3.) depending on the 72 !! coordinate chosen : 73 !! ln_zco=T z-coordinate (forced if lk_zco) 74 !! ln_zps=T z-coordinate with partial steps 75 !! ln_zco=T s-coordinate 76 !! 77 !! ** Action : define gdep., e3., mbathy and bathy 78 !!---------------------------------------------------------------------- 79 INTEGER :: ioptio = 0 ! temporary integer 80 !! 70 81 NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco 71 82 !!---------------------------------------------------------------------- 72 83 73 ! Read Namelist nam_zgr : vertical coordinate' 74 ! --------------------- 75 REWIND ( numnam ) 84 REWIND ( numnam ) ! Read Namelist nam_zgr : vertical coordinate' 76 85 READ ( numnam, nam_zgr ) 77 86 78 ! Parameter control and print 79 ! --------------------------- 80 ! Control print 81 IF(lwp) THEN 87 IF(lwp) THEN ! Control print 82 88 WRITE(numout,*) 83 89 WRITE(numout,*) 'dom_zgr : vertical coordinate' … … 89 95 ENDIF 90 96 91 ! Check Vertical coordinate options 92 ioptio = 0 97 ioptio = 0 ! Check Vertical coordinate options 93 98 IF( ln_zco ) ioptio = ioptio + 1 94 99 IF( ln_zps ) ioptio = ioptio + 1 95 100 IF( ln_sco ) ioptio = ioptio + 1 96 IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 97 101 IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 98 102 IF( lk_zco ) THEN 99 103 IF(lwp) WRITE(numout,*) ' z-coordinate with reduced incore memory requirement' 100 IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' )104 IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' ) 101 105 ENDIF 102 106 103 107 ! Build the vertical coordinate system 104 108 ! ------------------------------------ 105 106 109 CALL zgr_z ! Reference z-coordinate system (always called) 107 108 110 CALL zgr_bat ! Bathymetry fields (levels and meters) 109 110 111 IF( ln_zco ) CALL zgr_zco ! z-coordinate 111 112 112 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 113 114 113 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 115 116 !!bug gm control print: 117 IF( nprint == 1 .AND. lwp ) THEN 118 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 119 WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), & 120 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) ) 121 WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), & 122 & ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), & 123 & ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), & 124 & ' w ', MINVAL( fse3w(:,:,:) ) 125 126 WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), & 127 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) ) 128 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), & 129 & ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), & 130 & ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), & 131 & ' w ', MAXVAL( fse3w(:,:,:) ) 132 ENDIF 133 !!!bug gm 134 114 ! 135 115 END SUBROUTINE dom_zgr 136 116 … … 153 133 !! 154 134 !! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m) 155 !! - e3t_0, e3w_0 : scale factors at T- and W-levels (m) 156 !! 157 !! Reference : 158 !! Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 159 !! 160 !! History : 161 !! 9.0 ! 03-08 (G. Madec) F90: Free form and module 162 !!---------------------------------------------------------------------- 163 !! * Local declarations 135 !! - e3t_0 , e3w_0 : scale factors at T- and W-levels (m) 136 !! 137 !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 138 !!---------------------------------------------------------------------- 164 139 INTEGER :: jk ! dummy loop indices 165 140 REAL(wp) :: zt, zw ! temporary scalars 166 REAL(wp) :: & 167 zsur , za0, za1, zkth, zacr, & ! Values set from parameters in 168 zdzmin, zhmax ! par_CONFIG_Rxx.h90 141 REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in 142 REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90 169 143 !!---------------------------------------------------------------------- 170 144 … … 177 151 ! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 178 152 ! 179 IF(ppa1 == pp_to_be_computed .AND. &153 IF( ppa1 == pp_to_be_computed .AND. & 180 154 & ppa0 == pp_to_be_computed .AND. & 181 155 & ppsur == pp_to_be_computed ) THEN 182 za1 = ( ppdzmin - pphmax / FLOAT(jpk-1) ) & 183 / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) & 184 & * ( LOG( COSH( (jpk - ppkth) / ppacr) ) & 185 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 186 187 za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 188 zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) 189 190 ELSE 156 ! 157 za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) & 158 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) & 159 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 160 za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 161 zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) 162 ELSE 191 163 za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur 192 ENDIF 193 194 195 ! Parameter print 196 ! --------------- 197 IF(lwp) THEN 164 ENDIF 165 166 IF(lwp) THEN ! Parameter print 198 167 WRITE(numout,*) 199 168 WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates' … … 222 191 ! ====================== 223 192 IF( ppkth == 0.e0 ) THEN ! uniform vertical grid 224 225 193 za1 = zhmax / FLOAT(jpk-1) 226 194 DO jk = 1, jpk … … 232 200 e3t_0 (jk) = za1 233 201 END DO 234 235 ELSE 236 202 ELSE ! Madec & Imbard 1996 function 237 203 DO jk = 1, jpk 238 204 zw = FLOAT( jk ) 239 205 zt = FLOAT( jk ) + 0.5 240 gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) ) ) 241 gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) ) ) 242 e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth)/zacr ) 243 e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth)/zacr ) 244 END DO 245 gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero 246 247 ENDIF 248 249 ! Control and print 250 ! ================== 251 252 IF(lwp) THEN 206 gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) ) 207 gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) ) 208 e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth) / zacr ) 209 e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth) / zacr ) 210 END DO 211 gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero 212 ENDIF 213 214 IF(lwp) THEN ! control print 253 215 WRITE(numout,*) 254 216 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' … … 256 218 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk ) 257 219 ENDIF 258 259 DO jk = 1, jpk 260 IF( e3w_0(jk) <= 0. .OR. e3t_0(jk) <= 0. ) CALL ctl_stop( ' e3w or e3t =< 0 ' ) 261 IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0. ) CALL ctl_stop( ' gdepw or gdept < 0 ' ) 262 END DO 263 220 DO jk = 1, jpk ! control positivity 221 IF( e3w_0 (jk) <= 0.e0 .OR. e3t_0 (jk) <= 0.e0 ) CALL ctl_stop( ' e3w or e3t =< 0 ' ) 222 IF( gdepw_0(jk) < 0.e0 .OR. gdept_0(jk) < 0.e0 ) CALL ctl_stop( ' gdepw or gdept < 0 ' ) 223 END DO 224 ! 264 225 END SUBROUTINE zgr_z 265 226 … … 298 259 !! ** Action : - mbathy: level bathymetry (in level index) 299 260 !! - bathy : meter bathymetry (in meters) 300 !! 301 !! History : 302 !! 9.0 ! 03-08 (G. Madec) Original code 303 !! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates 304 !!---------------------------------------------------------------------- 305 !! * Modules used 261 !!---------------------------------------------------------------------- 306 262 USE iom 307 308 !! * Local declarations 309 INTEGER :: ji, jj, jl, jk ! dummy loop indices 310 INTEGER :: inum ! temporary logical unit 311 INTEGER :: & 312 ii_bump, ij_bump, ih ! bump center position 313 INTEGER , DIMENSION(jpidta,jpjdta) :: & 314 idta ! global domain integer data 315 REAL(wp) :: & 316 r_bump, h_bump, h_oce, & ! bump characteristics 317 zi, zj, zh ! temporary scalars 318 REAL(wp), DIMENSION(jpidta,jpjdta) :: & 319 zdta ! global domain scalar data 263 !! 264 INTEGER :: ji, jj, jl, jk ! dummy loop indices 265 INTEGER :: inum ! temporary logical unit 266 INTEGER :: ii_bump, ij_bump, ih ! bump center position 267 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 268 REAL(wp) :: zi , zj , zh ! temporary scalars 269 INTEGER , DIMENSION(jpidta,jpjdta) :: idta ! global domain integer data 270 REAL(wp), DIMENSION(jpidta,jpjdta) :: zdta ! global domain scalar data 320 271 !!---------------------------------------------------------------------- 321 272 … … 324 275 IF(lwp) WRITE(numout,*) ' ~~~~~~~' 325 276 326 ! ======================================== 327 ! global domain level and meter bathymetry (idta,zdta) 328 ! ======================================== 329 ! ! =============== ! 330 IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand ! 331 ! ! =============== ! 332 277 ! ! ================== ! 278 IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand ! 279 ! ! ================== ! 280 ! ! global domain level and meter bathymetry (idta,zdta) 281 ! 333 282 IF( ntopo == 0 ) THEN ! flat basin 334 335 283 IF(lwp) WRITE(numout,*) 336 284 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin' 337 338 idta(:,:) = jpkm1 ! flat basin 339 zdta(:,:) = gdepw_0(jpk) 285 idta(:,:) = jpkm1 ! before last level 286 zdta(:,:) = gdepw_0(jpk) ! last w-point depth 340 287 h_oce = gdepw_0(jpk) 341 342 ELSE ! bump 288 ELSE ! bump centered in the basin 343 289 IF(lwp) WRITE(numout,*) 344 290 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump' 345 346 ii_bump = jpidta / 2 ! i-index of the bump center 347 ij_bump = jpjdta / 2 ! j-index of the bump center 348 r_bump = 50000.e0 ! bump radius (meters) 349 h_bump = 2700.e0 ! bump height (meters) 350 h_oce = gdepw_0(jpk) ! background ocean depth (meters) 291 ii_bump = jpidta / 2 ! i-index of the bump center 292 ij_bump = jpjdta / 2 ! j-index of the bump center 293 r_bump = 50000.e0 ! bump radius (meters) 294 h_bump = 2700.e0 ! bump height (meters) 295 h_oce = gdepw_0(jpk) ! background ocean depth (meters) 351 296 IF(lwp) WRITE(numout,*) ' bump characteristics: ' 352 297 IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump … … 354 299 IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index' 355 300 IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters' 356 ! zdta :357 DO jj = 1, jpjdta 301 ! 302 DO jj = 1, jpjdta ! zdta : 358 303 DO ji = 1, jpidta 359 304 zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump … … 362 307 END DO 363 308 END DO 364 365 ! idta : 366 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 309 ! ! idta : 310 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 367 311 idta(:,:) = jpkm1 368 ELSE ! z-coordinate (zco or zps): step-like topography312 ELSE ! z-coordinate (zco or zps): step-like topography 369 313 idta(:,:) = jpkm1 370 314 DO jk = 1, jpkm1 … … 377 321 ENDIF 378 322 ENDIF 379 380 ! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio)323 ! ! set GLOBAL boundary conditions 324 ! ! Caution : idta on the global domain: use of jperio, not nperio 381 325 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 382 326 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1.e0 … … 396 340 ENDIF 397 341 398 ! ======================================= 399 ! local domain level and meter bathymetry (mbathy,bathy) 400 ! ======================================= 401 402 mbathy(:,:) = 0 ! set to zero extra halo points 403 bathy (:,:) = 0.e0 ! (require for mpp case) 404 405 DO jj = 1, nlcj ! interior values 342 ! ! local domain level and meter bathymetries (mbathy,bathy) 343 mbathy(:,:) = 0 ! set to zero extra halo points 344 bathy (:,:) = 0.e0 ! (require for mpp case) 345 DO jj = 1, nlcj ! interior values 406 346 DO ji = 1, nlci 407 347 mbathy(ji,jj) = idta( mig(ji), mjg(jj) ) … … 409 349 END DO 410 350 END DO 411 412 ! ! =============== !413 ELSEIF( ntopo == 1 ) THEN ! read in file !414 ! ! =============== !415 416 IF( ln_zco ) THEN 417 CALL iom_open ( 'bathy_level.nc', inum ) ! Level bathymetry351 ! 352 ! ! ================ ! 353 ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain) 354 ! ! ================ ! 355 ! 356 IF( ln_zco ) THEN ! zco : read level bathymetry 357 CALL iom_open( 'bathy_level.nc', inum ) 418 358 CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy ) 419 359 CALL iom_close (inum) 420 360 mbathy(:,:) = INT( bathy(:,:) ) 421 361 ENDIF 422 423 IF( ln_zps .OR. ln_sco ) THEN 424 CALL iom_open ( 'bathy_meter.nc', inum ) ! meter bathymetry 362 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 363 CALL iom_open( 'bathy_meter.nc', inum ) 425 364 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 426 365 CALL iom_close (inum) … … 429 368 ELSE ! error ! 430 369 ! ! =============== ! 431 WRITE(ctmp1,*) ' 370 WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo 432 371 CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) ) 433 372 ENDIF 434 435 ! ======================= 436 ! NO closed seas or lakes 437 ! ======================= 438 439 IF( nclosea == 0 ) THEN 440 DO jl = 1, jpncs 373 ! 374 ! ! =========================== ! 375 IF( nclosea == 0 ) THEN ! NO closed seas or lakes ! 376 DO jl = 1, jpncs ! =========================== ! 441 377 DO jj = ncsj1(jl), ncsj2(jl) 442 378 DO ji = ncsi1(jl), ncsi2(jl) 443 mbathy(ji,jj) = 0 ! suppress closed seas 444 bathy (ji,jj) = 0.e0 ! and lakes 445 !!bug gm bathy (ji,jj) = 300.e0 ! and lakes 379 mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry 380 bathy (ji,jj) = 0.e0 446 381 END DO 447 382 END DO 448 383 END DO 449 384 ENDIF 450 451 385 #if defined key_orca_lev10 452 ! 10 time the vertical resolution 453 mbathy(:,:) = 10 * mbathy(:,:) 454 IF(lwp) WRITE(numout,*) ' ATTENTION: 300 niveaux avec bathy levels "vraie?"' 386 ! ! ================= ! 387 mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 ! 388 ! ! ================= ! 389 IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' ) 455 390 #endif 456 ! =========== 457 ! Zoom domain 458 ! =========== 459 460 IF( lzoom ) CALL zgr_bat_zoom 461 462 ! ================ 463 ! Bathymetry check 464 ! ================ 465 466 CALL zgr_bat_ctl 467 391 ! ! =============== ! 392 IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain ! 393 ! ! =============== ! 394 395 ! ! =================== ! 396 IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check ! 397 ! ! =================== ! 468 398 END SUBROUTINE zgr_bat 469 399 … … 479 409 !! 480 410 !! ** Action : - update mbathy: level bathymetry (in level index) 481 !! 482 !! History : 483 !! 9.0 ! 03-08 (G. Madec) Original code 484 !!---------------------------------------------------------------------- 485 !! * Local variables 411 !!---------------------------------------------------------------------- 486 412 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 487 413 !!---------------------------------------------------------------------- 488 414 ! 489 415 IF(lwp) WRITE(numout,*) 490 416 IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain' 491 417 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~' 492 418 ! 493 419 ! Zoom domain 494 420 ! =========== 495 421 ! 496 422 ! Forced closed boundary if required 497 IF( lzoom_ w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : )= 0498 IF( lzoom_ s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom)) = 0499 IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0500 IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0501 423 IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0 424 IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0 425 IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0 426 IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0 427 ! 502 428 ! Configuration specific domain modifications 503 429 ! (here, ORCA arctic configuration: suppress Med Sea) … … 508 434 ! ! ======================= 509 435 IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea' 510 511 436 ii0 = 141 ; ii1 = 162 ! Sea box i,j indices 512 437 ij0 = 98 ; ij1 = 110 … … 515 440 ! ! ======================= 516 441 IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea' 517 518 442 ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe 519 443 ij0 = 314 ; ij1 = 370 520 521 444 END SELECT 522 445 ! … … 524 447 ! 525 448 ENDIF 526 449 ! 527 450 END SUBROUTINE zgr_bat_zoom 528 451 … … 549 472 !! ** Action : - update mbathy: level bathymetry (in level index) 550 473 !! - update bathy : meter bathymetry (in meters) 551 !! 552 !! History : 553 !! 9.0 ! 03-08 (G. Madec) Original code 554 !! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates 555 !!---------------------------------------------------------------------- 556 !! * Local declarations 557 INTEGER :: ji, jj, jl ! dummy loop indices 558 INTEGER :: & 559 icompt, ibtest, ikmax ! temporary integers 560 REAL(wp), DIMENSION(jpi,jpj) :: & 561 zbathy ! temporary workspace 474 !!---------------------------------------------------------------------- 475 INTEGER :: ji, jj, jl ! dummy loop indices 476 INTEGER :: icompt, ibtest, ikmax ! temporary integers 477 REAL(wp), DIMENSION(jpi,jpj) :: zbathy ! temporary workspace 562 478 !!---------------------------------------------------------------------- 563 479 … … 566 482 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 567 483 568 ! ================ 569 ! Bathymetry check 570 ! ================ 571 572 IF( .NOT. lk_c1d ) THEN 573 574 ! Suppress isolated ocean grid points 575 484 ! ! Suppress isolated ocean grid points 485 IF(lwp) WRITE(numout,*) 486 IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points' 487 IF(lwp) WRITE(numout,*)' -----------------------------------' 488 icompt = 0 489 DO jl = 1, 2 490 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 491 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 492 mbathy(jpi,:) = mbathy( 2 ,:) 493 ENDIF 494 DO jj = 2, jpjm1 495 DO ji = 2, jpim1 496 ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), & 497 & mbathy(ji,jj-1), mbathy(ji,jj+1) ) 498 IF( ibtest < mbathy(ji,jj) ) THEN 499 IF(lwp) WRITE(numout,*) ' the number of ocean level at ', & 500 & 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest 501 mbathy(ji,jj) = ibtest 502 icompt = icompt + 1 503 ENDIF 504 END DO 505 END DO 506 END DO 507 IF( icompt == 0 ) THEN 508 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' 509 ELSE 510 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' 511 ENDIF 512 IF( lk_mpp ) THEN 513 zbathy(:,:) = FLOAT( mbathy(:,:) ) 514 CALL lbc_lnk( zbathy, 'T', 1. ) 515 mbathy(:,:) = INT( zbathy(:,:) ) 516 ENDIF 517 518 ! ! East-west cyclic boundary conditions 519 IF( nperio == 0 ) THEN 520 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio 521 IF( lk_mpp ) THEN 522 IF( nbondi == -1 .OR. nbondi == 2 ) THEN 523 IF( jperio /= 1 ) mbathy(1,:) = 0 524 ENDIF 525 IF( nbondi == 1 .OR. nbondi == 2 ) THEN 526 IF( jperio /= 1 ) mbathy(nlci,:) = 0 527 ENDIF 528 ELSE 529 IF( ln_zco .OR. ln_zps ) THEN 530 mbathy( 1 ,:) = 0 531 mbathy(jpi,:) = 0 532 ELSE 533 mbathy( 1 ,:) = jpkm1 534 mbathy(jpi,:) = jpkm1 535 ENDIF 536 ENDIF 537 ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 538 IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio 539 mbathy( 1 ,:) = mbathy(jpim1,:) 540 mbathy(jpi,:) = mbathy( 2 ,:) 541 ELSEIF( nperio == 2 ) THEN 542 IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio 543 ELSE 544 IF(lwp) WRITE(numout,*) ' e r r o r' 545 IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio 546 ! STOP 'dom_mba' 547 ENDIF 548 549 ! Set to zero mbathy over islands if necessary (lk_isl=F) 550 IF( .NOT. lk_isl ) THEN ! No island 576 551 IF(lwp) WRITE(numout,*) 577 IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points' 578 IF(lwp) WRITE(numout,*)' -----------------------------------' 579 580 icompt = 0 581 582 DO jl = 1, 2 583 584 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 585 mbathy( 1 ,:) = mbathy(jpim1,:) 586 mbathy(jpi,:) = mbathy( 2 ,:) 587 ENDIF 588 DO jj = 2, jpjm1 589 DO ji = 2, jpim1 590 ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), & 591 & mbathy(ji,jj-1),mbathy(ji,jj+1) ) 592 IF( ibtest < mbathy(ji,jj) ) THEN 593 IF(lwp) WRITE(numout,*) ' the number of ocean level at ', & 594 & 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest 595 mbathy(ji,jj) = ibtest 596 icompt = icompt + 1 597 ENDIF 598 END DO 599 END DO 600 601 END DO 602 IF( icompt == 0 ) THEN 603 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' 604 ELSE 605 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' 606 ENDIF 607 IF( lk_mpp ) THEN 552 IF(lwp) WRITE(numout,*) ' mbathy set to 0 over islands' 553 IF(lwp) WRITE(numout,*) ' ----------------------------' 554 ! 555 mbathy(:,:) = MAX( 0, mbathy(:,:) ) 556 ! 557 ! Boundary condition on mbathy 558 IF( .NOT.lk_mpp ) THEN 559 !!gm !!bug ??? think about it ! 560 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab 608 561 zbathy(:,:) = FLOAT( mbathy(:,:) ) 609 562 CALL lbc_lnk( zbathy, 'T', 1. ) 610 563 mbathy(:,:) = INT( zbathy(:,:) ) 611 564 ENDIF 612 613 ! 3.2 East-west cyclic boundary conditions614 615 IF( nperio == 0 ) THEN616 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west', &617 ' boundary: nperio = ', nperio618 IF( lk_mpp ) THEN619 IF( nbondi == -1 .OR. nbondi == 2 ) THEN620 IF( jperio /= 1 ) mbathy(1,:) = 0621 ENDIF622 IF( nbondi == 1 .OR. nbondi == 2 ) THEN623 IF( jperio /= 1 ) mbathy(nlci,:) = 0624 ENDIF625 ELSE626 IF( ln_zco .OR. ln_zps ) THEN627 mbathy( 1 ,:) = 0628 mbathy(jpi,:) = 0629 ELSE630 mbathy( 1 ,:) = jpkm1631 mbathy(jpi,:) = jpkm1632 ENDIF633 ENDIF634 ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN635 IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions', &636 ' on mbathy: nperio = ', nperio637 mbathy( 1 ,:) = mbathy(jpim1,:)638 mbathy(jpi,:) = mbathy( 2 ,:)639 ELSEIF( nperio == 2 ) THEN640 IF(lwp) WRITE(numout,*) ' equatorial boundary conditions', &641 ' on mbathy: nperio = ', nperio642 ELSE643 IF(lwp) WRITE(numout,*) ' e r r o r'644 IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio645 ! STOP 'dom_mba'646 ENDIF647 648 ! Set to zero mbathy over islands if necessary (lk_isl=F)649 IF( .NOT. lk_isl ) THEN ! No island650 IF(lwp) WRITE(numout,*)651 IF(lwp) WRITE(numout,*) ' mbathy set to 0 over islands'652 IF(lwp) WRITE(numout,*) ' ----------------------------'653 654 mbathy(:,:) = MAX( 0, mbathy(:,:) )655 656 ! Boundary condition on mbathy657 IF( .NOT.lk_mpp ) THEN658 !!bug ??? y reflechir!659 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab660 zbathy(:,:) = FLOAT( mbathy(:,:) )661 CALL lbc_lnk( zbathy, 'T', 1. )662 mbathy(:,:) = INT( zbathy(:,:) )663 ENDIF664 665 ENDIF666 667 565 ENDIF 668 566 669 567 ! Number of ocean level inferior or equal to jpkm1 670 671 568 ikmax = 0 672 569 DO jj = 1, jpj … … 675 572 END DO 676 573 END DO 677 !!! test a faire: ikmax = MAX( mbathy(:,:) ) ??? 678 574 !!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ??? 679 575 IF( ikmax > jpkm1 ) THEN 680 576 IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1' … … 685 581 ENDIF 686 582 583 ! control print 687 584 IF( lwp .AND. nprint == 1 ) THEN 688 585 WRITE(numout,*) 689 WRITE(numout,*) ' bathymetric field '586 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels ' 690 587 WRITE(numout,*) ' ------------------' 691 WRITE(numout,*) ' number of non-zero T-levels ' 692 CALL prihin( mbathy, jpi, jpj, 1, jpi, & 693 1 , 1 , jpj, 1, 3 , & 694 numout ) 695 WRITE(numout,*) 696 ENDIF 697 588 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout ) 589 WRITE(numout,*) 590 ENDIF 591 ! 698 592 END SUBROUTINE zgr_bat_ctl 699 593 … … 706 600 !! 707 601 !! ** Method : set 3D coord. arrays to reference 1D array if lk_zco=T 708 !!709 !! History :710 !! ! 06-04 (R. Benshila, G. Madec) Original code711 602 !!---------------------------------------------------------------------- 712 603 INTEGER :: jk 713 604 !!---------------------------------------------------------------------- 714 605 ! 715 606 IF( .NOT.lk_zco ) THEN 716 607 DO jk = 1, jpk … … 727 618 END DO 728 619 ENDIF 729 620 ! 730 621 END SUBROUTINE zgr_zco 731 622 … … 787 678 !! - - - - - - - gdept, gdepw and e3. are positives 788 679 !! 789 !! Reference : 790 !! Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 791 !! 792 !! History : 793 !! 8.5 ! 02-09 (A. Bozec, G. Madec) F90: Free form and module 794 !! 9.0 ! 02-09 (A. de Miranda) rigid-lid + islands 795 !!---------------------------------------------------------------------- 796 !! * Local declarations 797 INTEGER :: ji, jj, jk ! dummy loop indices 798 INTEGER :: & 799 ik, it ! temporary integers 800 801 REAL(wp) :: & 802 ze3tp, ze3wp, & ! Last ocean level thickness at T- and W-points 803 zdepwp, & ! Ajusted ocean depth to avoid too small e3t 804 zdepth, & ! " " 805 zmax, zmin, & ! Maximum and minimum depth 806 zdiff ! temporary scalar 807 808 REAL(wp), DIMENSION(jpi,jpj) :: & 809 zprt ! " " 810 811 LOGICAL :: ll_print ! Allow control print for debugging 812 680 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 681 !!---------------------------------------------------------------------- 682 INTEGER :: ji, jj, jk ! dummy loop indices 683 INTEGER :: ik, it ! temporary integers 684 LOGICAL :: ll_print ! Allow control print for debugging 685 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 686 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 687 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 688 REAL(wp) :: zdiff ! temporary scalar 689 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprt ! 3D workspace 813 690 !!--------------------------------------------------------------------- 814 !! OPA8.5, LODYC-IPSL (2002)815 !!---------------------------------------------------------------------816 817 ! Local variable for debugging818 ll_print=.FALSE.819 !!! ll_print=.TRUE.820 821 ! Initialization of constant822 zmax = gdepw_0(jpk) + e3t_0(jpk)823 zmin = gdepw_0(4)824 825 ! Ocean depth826 IF(lwp .AND. ll_print) THEN827 WRITE(numout,*)828 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'829 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )830 ENDIF831 691 832 692 IF(lwp) WRITE(numout,*) … … 835 695 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 836 696 697 ll_print=.FALSE. ! Local variable for debugging 698 !! ll_print=.TRUE. 699 700 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth 701 WRITE(numout,*) 702 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 703 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 704 ENDIF 705 837 706 838 707 ! bathymetry in level (from bathy_meter) 839 708 ! =================== 840 841 ! initialize mbathy to the maximum ocean level available 842 mbathy(:,:) = jpkm1 843 844 ! storage of land and island's number (zera and negative values) in mbathy 709 zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) 710 zmin = gdepw_0(4) ! minimum depth = 3 levels 711 712 mbathy(:,:) = jpkm1 ! initialize mbathy to the maximum ocean level available 713 714 ! ! storage of land and island's number (zera and negative values) in mbathy 715 WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) ) 716 717 ! bounded value of bathy 718 !!gm bathy(:,:) = MIN( zmax, MAX( bathy(:,:), zmin ) ) !!gm : simpler as zmin is > 0 845 719 DO jj = 1, jpj 846 720 DO ji= 1, jpi 847 IF( bathy(ji,jj) <= 0. ) mbathy(ji,jj) = INT( bathy(ji,jj) ) 848 END DO 849 END DO 850 851 ! bounded value of bathy 852 ! minimum depth == 3 levels 853 ! maximum depth == gdepw_0(jpk)+e3t_0(jpk) 854 ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk) 855 DO jj = 1, jpj 856 DO ji= 1, jpi 857 IF( bathy(ji,jj) <= 0. ) THEN 858 bathy(ji,jj) = 0.e0 859 ELSE 860 bathy(ji,jj) = MAX( bathy(ji,jj), zmin ) 861 bathy(ji,jj) = MIN( bathy(ji,jj), zmax ) 721 IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0 722 ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) ) 862 723 ENDIF 863 724 END DO … … 877 738 END DO 878 739 879 ! Scale factors and depth at T- and W-points 880 881 ! intitialization to the reference z-coordinate 882 DO jk = 1, jpk 883 gdept(:,:,jk) = gdept_0(jk) 884 gdepw(:,:,jk) = gdepw_0(jk) 885 e3t (:,:,jk) = e3t_0 (jk) 886 e3w (:,:,jk) = e3w_0 (jk) 887 END DO 888 hdept(:,:) = gdept(:,:,2 ) 889 hdepw(:,:) = gdepw(:,:,3 ) 890 891 ! 892 DO jj = 1, jpj 893 DO ji = 1, jpi 894 ik = mbathy(ji,jj) 895 ! ocean point only 896 IF( ik > 0 ) THEN 897 ! max ocean level case 898 IF( ik == jpkm1 ) THEN 899 zdepwp = bathy(ji,jj) 900 ze3tp = bathy(ji,jj) - gdepw_0(ik) 901 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) ) 902 e3t(ji,jj,ik ) = ze3tp 903 e3t(ji,jj,ik+1) = ze3tp 904 e3w(ji,jj,ik ) = ze3wp 905 e3w(ji,jj,ik+1) = ze3tp 906 gdepw(ji,jj,ik+1) = zdepwp 907 gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp 908 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp 909 ! standard case 910 ELSE 911 !!alex 912 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN 913 gdepw(ji,jj,ik+1) = bathy(ji,jj) 914 ELSE 915 !!alex ctl write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1) 916 gdepw(ji,jj,ik+1) = gdepw_0(ik+1) 917 ENDIF 918 !!Alex 919 !!Alex gdepw(ji,jj,ik+1) = bathy(ji,jj) 920 !!bug gm verifier les gdepw_0 921 gdept(ji,jj,ik ) = gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) ) & 922 & * ((gdept_0( ik ) - gdepw_0(ik) ) & 923 & / ( gdepw_0( ik+1) - gdepw_0(ik) )) 924 e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 925 & / ( gdepw_0( ik+1) - gdepw_0(ik) ) 926 e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) & 927 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 928 ! ... on ik+1 929 e3w(ji,jj,ik+1) = e3t(ji,jj,ik) 930 e3t(ji,jj,ik+1) = e3t(ji,jj,ik) 931 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) 932 ENDIF 933 ENDIF 934 END DO 935 END DO 936 937 it = 0 938 DO jj = 1, jpj 939 DO ji = 1, jpi 940 ik = mbathy(ji,jj) 941 ! ocean point only 942 IF( ik > 0 ) THEN 943 ! bathymetry output 944 hdept(ji,jj) = gdept(ji,jj,ik ) 945 hdepw(ji,jj) = gdepw(ji,jj,ik+1) 946 e3tp (ji,jj) = e3t(ji,jj,ik ) 947 e3wp (ji,jj) = e3w(ji,jj,ik ) 948 ! test 949 zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) 950 IF( zdiff <= 0. .AND. lwp ) THEN 951 it=it+1 952 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 953 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 954 WRITE(numout,*) ' gdept= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1), & 955 ' zdiff = ', zdiff 956 WRITE(numout,*) ' e3tp = ', e3t(ji,jj,ik ), ' e3wp = ', e3w(ji,jj,ik ) 957 ENDIF 958 ENDIF 959 END DO 740 ! Scale factors and depth at T- and W-points 741 DO jk = 1, jpk ! intitialization to the reference z-coordinate 742 gdept(:,:,jk) = gdept_0(jk) 743 gdepw(:,:,jk) = gdepw_0(jk) 744 e3t (:,:,jk) = e3t_0 (jk) 745 e3w (:,:,jk) = e3w_0 (jk) 746 END DO 747 hdept(:,:) = gdept(:,:,2 ) 748 hdepw(:,:) = gdepw(:,:,3 ) 749 ! 750 DO jj = 1, jpj 751 DO ji = 1, jpi 752 ik = mbathy(ji,jj) 753 IF( ik > 0 ) THEN ! ocean point only 754 ! max ocean level case 755 IF( ik == jpkm1 ) THEN 756 zdepwp = bathy(ji,jj) 757 ze3tp = bathy(ji,jj) - gdepw_0(ik) 758 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) ) 759 e3t(ji,jj,ik ) = ze3tp 760 e3t(ji,jj,ik+1) = ze3tp 761 e3w(ji,jj,ik ) = ze3wp 762 e3w(ji,jj,ik+1) = ze3tp 763 gdepw(ji,jj,ik+1) = zdepwp 764 gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp 765 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp 766 ! 767 ELSE ! standard case 768 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj) 769 ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1) 770 ENDIF 771 !gm Bug? check the gdepw_0 772 ! ... on ik 773 gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 774 & * ((gdept_0( ik ) - gdepw_0(ik) ) & 775 & / ( gdepw_0( ik+1) - gdepw_0(ik) )) 776 e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 777 & / ( gdepw_0( ik+1) - gdepw_0(ik) ) 778 e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) & 779 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 780 ! ... on ik+1 781 e3w (ji,jj,ik+1) = e3t (ji,jj,ik) 782 e3t (ji,jj,ik+1) = e3t (ji,jj,ik) 783 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) 784 ENDIF 785 ENDIF 786 END DO 787 END DO 788 ! 789 it = 0 790 DO jj = 1, jpj 791 DO ji = 1, jpi 792 ik = mbathy(ji,jj) 793 IF( ik > 0 ) THEN ! ocean point only 794 hdept(ji,jj) = gdept(ji,jj,ik ) 795 hdepw(ji,jj) = gdepw(ji,jj,ik+1) 796 e3tp (ji,jj) = e3t(ji,jj,ik ) 797 e3wp (ji,jj) = e3w(ji,jj,ik ) 798 ! test 799 zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) 800 IF( zdiff <= 0. .AND. lwp ) THEN 801 it = it + 1 802 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 803 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 804 WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff 805 WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik ) 806 ENDIF 807 ENDIF 808 END DO 960 809 END DO 961 810 962 811 ! Scale factors and depth at U-, V-, UW and VW-points 963 964 ! initialisation to z-scale factors 965 DO jk = 1, jpk 812 DO jk = 1, jpk ! initialisation to z-scale factors 966 813 e3u (:,:,jk) = e3t_0(jk) 967 814 e3v (:,:,jk) = e3t_0(jk) … … 969 816 e3vw(:,:,jk) = e3w_0(jk) 970 817 END DO 971 972 ! Computed as the minimum of neighbooring scale factors 973 DO jk = 1,jpk 974 DO jj = 1, jpjm1 975 DO ji = 1, fs_jpim1 ! vector opt. 976 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 977 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk)) 978 e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 979 e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 980 END DO 981 END DO 982 END DO 818 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 819 DO jj = 1, jpjm1 820 DO ji = 1, fs_jpim1 ! vector opt. 821 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 822 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk)) 823 e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 824 e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 825 END DO 826 END DO 827 END DO 828 ! ! Boundary conditions 829 CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. ) 830 CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. ) 831 ! 832 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 833 WHERE( e3u (:,:,jk) == 0.e0 ) e3u (:,:,jk) = e3t_0(jk) 834 WHERE( e3v (:,:,jk) == 0.e0 ) e3v (:,:,jk) = e3t_0(jk) 835 WHERE( e3uw(:,:,jk) == 0.e0 ) e3uw(:,:,jk) = e3w_0(jk) 836 WHERE( e3vw(:,:,jk) == 0.e0 ) e3vw(:,:,jk) = e3w_0(jk) 837 END DO 838 839 ! Scale factor at F-point 840 DO jk = 1, jpk ! initialisation to z-scale factors 841 e3f (:,:,jk) = e3t_0(jk) 842 END DO 843 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 844 DO jj = 1, jpjm1 845 DO ji = 1, fs_jpim1 ! vector opt. 846 e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) 847 END DO 848 END DO 849 END DO 850 CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions 851 ! 852 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 853 WHERE( e3f(:,:,jk) == 0.e0 ) e3f(:,:,jk) = e3t_0(jk) 854 END DO 855 !!gm bug ? : must be a do loop with mj0,mj1 856 ! 857 e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 858 e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 859 e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 860 e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 861 e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 862 863 ! Control of the sign 864 IF( MINVAL( e3t (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' ) 865 IF( MINVAL( e3w (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' ) 866 IF( MINVAL( gdept(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) 867 IF( MINVAL( gdepw(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) 983 868 984 ! Boundary conditions 985 CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. ) 986 CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. ) 987 988 ! set to z-scale factor if zero (i.e. along closed boundaries) 989 DO jk = 1, jpk 990 DO jj = 1, jpj 991 DO ji = 1, jpi 992 IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk) 993 IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk) 994 IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk) 995 IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk) 996 END DO 997 END DO 998 END DO 999 1000 ! Scale factor at F-point 1001 1002 ! initialisation to z-scale factors 1003 DO jk = 1, jpk 1004 e3f (:,:,jk) = e3t_0(jk) 1005 END DO 1006 1007 ! Computed as the minimum of neighbooring V-scale factors 1008 DO jk = 1, jpk 1009 DO jj = 1, jpjm1 1010 DO ji = 1, fs_jpim1 ! vector opt. 1011 e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) 1012 END DO 1013 END DO 1014 END DO 1015 ! Boundary conditions 1016 CALL lbc_lnk( e3f, 'F', 1. ) 1017 1018 ! set to z-scale factor if zero (i.e. along closed boundaries) 1019 DO jk = 1, jpk 1020 DO jj = 1, jpj 1021 DO ji = 1, jpi 1022 IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk) 1023 END DO 1024 END DO 1025 END DO 1026 !!bug ? gm: must be a do loop with mj0,mj1 1027 1028 ! we duplicate factor scales for jj = 1 and jj = 2 1029 e3t(:,mj0(1),:) = e3t(:,mj0(2),:) 1030 e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 1031 e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 1032 e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 1033 e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 1034 1035 1036 1037 ! Compute gdep3w (vertical sum of e3w) 1038 1039 gdep3w (:,:,1) = 0.5 * e3w (:,:,1) 1040 DO jk = 2, jpk 1041 gdep3w (:,:,jk) = gdep3w (:,:,jk-1) + e3w (:,:,jk) 1042 END DO 869 ! Compute gdep3w (vertical sum of e3w) 870 gdep3w(:,:,1) = 0.5 * e3w(:,:,1) 871 DO jk = 2, jpk 872 gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 873 END DO 1043 874 1044 ! Control print 1045 9600 FORMAT(9x,' level gdept gdepw e3t e3w ') 1046 9610 FORMAT(10x,i4,4f9.2) 1047 IF(lwp .AND. ll_print) THEN 875 ! ! ================= ! 876 IF(lwp .AND. ll_print) THEN ! Control print ! 877 ! ! ================= ! 1048 878 DO jj = 1,jpj 1049 879 DO ji = 1, jpi 1050 ik = MAX(mbathy(ji,jj),1) 1051 zprt(ji,jj) = e3t(ji,jj,ik) 1052 END DO 1053 END DO 1054 WRITE(numout,*) 1055 WRITE(numout,*) 'domzgr e3t(mbathy)' 1056 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1057 DO jj = 1,jpj 1058 DO ji = 1, jpi 1059 ik = MAX(mbathy(ji,jj),1) 1060 zprt(ji,jj) = e3w(ji,jj,ik) 1061 END DO 1062 END DO 1063 WRITE(numout,*) 1064 WRITE(numout,*) 'domzgr e3w(mbathy)' 1065 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1066 DO jj = 1,jpj 1067 DO ji = 1, jpi 1068 ik = MAX(mbathy(ji,jj),1) 1069 zprt(ji,jj) = e3u(ji,jj,ik) 1070 END DO 1071 END DO 1072 1073 WRITE(numout,*) 1074 WRITE(numout,*) 'domzgr e3u(mbathy)' 1075 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1076 DO jj = 1,jpj 1077 DO ji = 1, jpi 1078 ik = MAX(mbathy(ji,jj),1) 1079 zprt(ji,jj) = e3v(ji,jj,ik) 1080 END DO 1081 END DO 1082 WRITE(numout,*) 1083 WRITE(numout,*) 'domzgr e3v(mbathy)' 1084 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1085 DO jj = 1,jpj 1086 DO ji = 1, jpi 1087 ik = MAX(mbathy(ji,jj),1) 1088 zprt(ji,jj) = e3f(ji,jj,ik) 1089 END DO 1090 END DO 1091 1092 WRITE(numout,*) 1093 WRITE(numout,*) 'domzgr e3f(mbathy)' 1094 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1095 DO jj = 1,jpj 1096 DO ji = 1, jpi 1097 ik = MAX(mbathy(ji,jj),1) 1098 zprt(ji,jj) = gdep3w(ji,jj,ik) 1099 END DO 1100 END DO 1101 WRITE(numout,*) 1102 WRITE(numout,*) 'domzgr gdep3w(mbathy)' 1103 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1104 880 ik = MAX( mbathy(ji,jj), 1 ) 881 zprt(ji,jj,1) = e3t (ji,jj,ik) 882 zprt(ji,jj,2) = e3w (ji,jj,ik) 883 zprt(ji,jj,3) = e3u (ji,jj,ik) 884 zprt(ji,jj,4) = e3v (ji,jj,ik) 885 zprt(ji,jj,5) = e3f (ji,jj,ik) 886 zprt(ji,jj,6) = gdep3w(ji,jj,ik) 887 END DO 888 END DO 889 WRITE(numout,*) 890 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 891 WRITE(numout,*) 892 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 893 WRITE(numout,*) 894 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 895 WRITE(numout,*) 896 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 897 WRITE(numout,*) 898 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 899 WRITE(numout,*) 900 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1105 901 ENDIF 1106 1107 1108 DO jk = 1,jpk 1109 DO jj = 1, jpj 1110 DO ji = 1, jpi 1111 IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN 1112 IF(lwp) THEN 1113 WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 ' 1114 WRITE(numout,*) ' ========= --------------- ' 1115 WRITE(numout,*) 1116 ENDIF 1117 STOP 'domzgr.psteps' 1118 ENDIF 1119 IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN 1120 IF(lwp) THEN 1121 WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 ' 1122 WRITE(numout,*) ' ========= ------------------ ' 1123 WRITE(numout,*) 1124 ENDIF 1125 STOP 'domzgr.psteps' 1126 ENDIF 1127 END DO 1128 END DO 1129 END DO 1130 1131 IF(lwp) THEN 1132 WRITE(numout,*) ' e3t lev 21 ' 1133 CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1134 WRITE(numout,*) ' e3w lev 21 ' 1135 CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1136 WRITE(numout,*) ' e3u lev 21 ' 1137 CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1138 WRITE(numout,*) ' e3v lev 21 ' 1139 CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1140 WRITE(numout,*) ' e3f lev 21 ' 1141 CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1142 WRITE(numout,*) ' e3t lev 22 ' 1143 CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1144 WRITE(numout,*) ' e3w lev 22 ' 1145 CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1146 WRITE(numout,*) ' e3u lev 22 ' 1147 CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1148 WRITE(numout,*) ' e3v lev 22 ' 1149 CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1150 WRITE(numout,*) ' e3f lev 22 ' 1151 CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1152 ENDIF 1153 1154 ! =========== 1155 ! Zoom domain 1156 ! =========== 1157 1158 IF( lzoom ) CALL zgr_bat_zoom 1159 1160 ! ================ 1161 ! Bathymetry check 1162 ! ================ 1163 1164 CALL zgr_bat_ctl 1165 902 903 ! ! =============== ! 904 IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain ! 905 ! ! =============== ! 906 907 ! ! =================== ! 908 IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check ! 909 ! ! =================== ! 1166 910 END SUBROUTINE zgr_zps 1167 911 … … 1223 967 !! schemes. 1224 968 !! 1225 !! Reference : 1226 !! Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 1227 !! 1228 !! History : 1229 !! ! 95-12 (G. Madec) Original code : s vertical coordinate 1230 !! ! 97-07 (G. Madec) lbc_lnk call 1231 !! ! 97-04 (J.-O. Beismann) 1232 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 1233 !! 9.0 ! 05-10 (A. Beckmann) new stretching function 1234 !!---------------------------------------------------------------------- 1235 !! * Local declarations 1236 INTEGER :: ji, jj, jk, jl 1237 INTEGER :: iip1, ijp1, iim1, ijm1 1238 REAL(wp) :: zrmax, ztaper 1239 1240 REAL(wp), DIMENSION(jpi,jpj) :: & 1241 zenv, ztmp, zmsk, zri, zrj, zhbat 1242 969 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 970 !!---------------------------------------------------------------------- 971 INTEGER :: ji, jj, jk, jl ! dummy loop argument 972 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 973 REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars 974 REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace 975 REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - - 976 !! 1243 977 NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max 1244 978 !!---------------------------------------------------------------------- 1245 979 1246 ! Read Namelist nam_zgr_sco : sigma-stretching parameters 1247 ! ------------------------- 1248 REWIND( numnam ) 980 REWIND( numnam ) ! Read Namelist nam_zgr_sco : sigma-stretching parameters 1249 981 READ ( numnam, nam_zgr_sco ) 1250 982 1251 IF(lwp) THEN 983 IF(lwp) THEN ! control print 1252 984 WRITE(numout,*) 1253 985 WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate' … … 1262 994 ENDIF 1263 995 1264 ! ??? 1265 ! ------------------ 1266 hift(:,:) = sbot_min 996 hift(:,:) = sbot_min ! set the minimum depth for the s-coordinate 1267 997 hifu(:,:) = sbot_min 1268 998 hifv(:,:) = sbot_min 1269 999 hiff(:,:) = sbot_min 1270 1271 1272 ! set maximum ocean depth 1273 ! ----------------------- 1274 DO jj = 1, jpj 1275 DO ji = 1, jpi 1276 bathy(ji,jj) = MIN( sbot_max, bathy(ji,jj) ) 1277 END DO 1278 END DO 1279 1280 ! Define the envelop bathymetry 1281 ! ============================= 1000 ! ! set maximum ocean depth 1001 bathy(:,:) = MIN( sbot_max, bathy(:,:) ) 1002 1003 ! ! ============================= 1004 ! ! Define the envelop bathymetry (hbatt) 1005 ! ! ============================= 1282 1006 ! Smooth the bathymetry (if required) 1283 1284 1007 scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea) 1285 1008 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1286 1287 1009 ! 1288 1010 ! use r-value to create hybrid coordinates 1289 1290 1011 DO jj = 1, jpj 1291 1012 DO ji = 1, jpi … … 1293 1014 END DO 1294 1015 END DO 1295 1296 1016 jl = 0 1297 1017 zrmax = 1.e0 … … 1300 1020 ! ! ================ ! 1301 1021 jl = jl + 1 1302 1303 1022 zrmax = 0.e0 1304 1023 zmsk(:,:) = 0.e0 1305 1306 1024 DO jj = 1, nlcj 1307 1025 DO ji = 1, nlci … … 1319 1037 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1320 1038 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 1321 ztmp(:,:) = zmsk(:,:) 1322 CALL lbc_lnk( zmsk, 'T', 1. ) 1039 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. ) 1323 1040 DO jj = 1, nlcj 1324 1041 DO ji = 1, nlci … … 1326 1043 END DO 1327 1044 END DO 1328 1329 WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )1330 1045 ! 1046 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 1047 ! 1331 1048 DO jj = 1, nlcj 1332 1049 DO ji = 1, nlci … … 1348 1065 END DO 1349 1066 END DO 1350 1067 ! 1351 1068 DO jj = 1, nlcj 1352 1069 DO ji = 1, nlci … … 1354 1071 END DO 1355 1072 END DO 1356 1073 ! 1357 1074 ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero 1358 ztmp(:,:) = zenv(:,:) 1359 CALL lbc_lnk( zenv, 'T', 1. ) 1075 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. ) 1360 1076 DO jj = 1, nlcj 1361 1077 DO ji = 1, nlci … … 1366 1082 END DO ! End loop ! 1367 1083 ! ! ================ ! 1368 1369 ! save envelop bathymetryin hbatt1084 ! 1085 ! ! envelop bathymetry saved in hbatt 1370 1086 hbatt(:,:) = zenv(:,:) 1371 1372 !!bug gm : CAUTION: tapering hard coded !!!! orca2 only 1373 !!bug gm NOT valid in mpp ===> taper have been changed 1374 1375 DO jj = 1, jpj1376 DO ji = 1, jpi1377 !!bug mpp ztaper = EXP( -(FLOAT(jj-74)/10.)**2 ) 1378 ztaper = EXP( -(gphit(ji,jj)/8)**2 )1379 hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)1380 END DO1381 END DO1382 1383 DO jj = 1, jpj1384 zrmax = EXP( -(gphit(10,jj)/8)**2)1385 ztaper = EXP( -(FLOAT(jj-74)/10.)**2)1386 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) jj, (FLOAT(jj-74)/10.), ztaper,(gphit(10,jj)/8), zrmax1387 END DO1388 1389 IF( nprint == 1 .AND. lwp ) THEN1390 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1391 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) 1392 ENDIF1393 1394 ! Control print1087 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN 1088 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 1089 DO jj = 1, jpj 1090 DO ji = 1, jpi 1091 ztaper = EXP( -(gphit(ji,jj)/8)**2 ) 1092 hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper) 1093 END DO 1094 END DO 1095 ENDIF 1096 ! 1097 IF(lwp) THEN ! Control print 1098 WRITE(numout,*) 1099 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 1100 WRITE(numout,*) 1101 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout ) 1102 IF( nprint == 1 ) THEN 1103 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) 1104 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) 1105 ENDIF 1106 ENDIF 1107 1108 ! ! ============================== 1109 ! ! hbatu, hbatv, hbatf fields 1110 ! ! ============================== 1395 1111 IF(lwp) THEN 1396 1112 WRITE(numout,*) 1397 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 1398 WRITE(numout,*) 1399 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout ) 1400 ENDIF 1401 1402 ! 1. computation of hbatu, hbatv, hbatf fields 1403 ! -------------------------------------------- 1404 1113 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min 1114 ENDIF 1405 1115 hbatu(:,:) = sbot_min 1406 1116 hbatv(:,:) = sbot_min 1407 1117 hbatf(:,:) = sbot_min 1408 1409 IF(lwp) THEN1410 WRITE(numout,*)1411 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min1412 WRITE(numout,*)1413 ENDIF1414 1415 1118 DO jj = 1, jpjm1 1416 1119 DO ji = 1, jpim1 1417 hbatu(ji,jj) = 0.5 *( hbatt(ji ,jj)+hbatt(ji+1,jj ) )1418 hbatv(ji,jj) = 0.5 *( hbatt(ji ,jj)+hbatt(ji ,jj+1) )1419 hbatf(ji,jj) = 0.25*( hbatt(ji ,jj)+hbatt(ji ,jj+1) &1420 +hbatt(ji+1,jj)+hbatt(ji+1,jj+1) )1120 hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 1121 hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 1122 hbatf(ji,jj) = 0.25* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 1123 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 1421 1124 END DO 1422 1125 END DO 1423 1126 ! 1424 1127 ! Apply lateral boundary condition 1425 ! CAUTION: retain non zero value in the initial file 1426 ! this should be OK for orca cfg, not for EEL 1427 zhbat(:,:) = hbatu(:,:) 1428 CALL lbc_lnk( hbatu, 'U', 1. ) 1128 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 1129 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. ) 1429 1130 DO jj = 1, jpj 1430 1131 DO ji = 1, jpi … … 1435 1136 END DO 1436 1137 END DO 1437 1438 zhbat(:,:) = hbatv(:,:) 1439 CALL lbc_lnk( hbatv, 'V', 1. ) 1138 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. ) 1440 1139 DO jj = 1, jpj 1441 1140 DO ji = 1, jpi … … 1446 1145 END DO 1447 1146 END DO 1448 1449 zhbat(:,:) = hbatf(:,:) 1450 CALL lbc_lnk( hbatf, 'F', 1. ) 1147 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. ) 1451 1148 DO jj = 1, jpj 1452 1149 DO ji = 1, jpi … … 1458 1155 END DO 1459 1156 1460 1461 1157 !!bug: key_helsinki a verifer 1462 1158 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) … … 1466 1162 1467 1163 IF( nprint == 1 .AND. lwp ) THEN 1468 WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff(:,:) ), &1469 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv(:,:) )1470 WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff(:,:) ), &1471 & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv(:,:) )1164 WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & 1165 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 1166 WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & 1167 & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 1472 1168 WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & 1473 1169 & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) … … 1477 1173 !! helsinki 1478 1174 1479 ! 2. Computation of gsig and esig fields1480 ! --------------------------------------1481 1482 ! 2.1 Coefficients for model level depth at w- and t-levels1483 1175 ! ! ======================= 1176 ! ! s-ccordinate fields (gdep., e3.) 1177 ! ! ======================= 1178 ! 1179 ! non-dimensional "sigma" for model level depth at w- and t-levels 1484 1180 DO jk = 1, jpk 1485 1181 gsigw(jk) = -fssig( FLOAT(jk)-0.5 ) 1486 1182 gsigt(jk) = -fssig( FLOAT(jk) ) 1487 1183 END DO 1488 1489 1184 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1490 1491 ! 2.2 Coefficients for vertical scale factors at w-, t- levels 1492 1493 !!bug gm : define it from analytical function, not like juste bellow.... 1494 !! or betteroffer the 2 possibilities.... 1495 DO jk=2,jpk 1496 esigw(jk)=gsigt(jk)-gsigt(jk-1) 1497 END DO 1498 DO jk=1,jpkm1 1499 esigt(jk)=gsigw(jk+1)-gsigw(jk) 1500 END DO 1501 esigw(1)=esigw(2) 1502 esigt(jpk)=esigt(jpkm1) 1503 1504 ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors 1505 1185 ! 1186 ! Coefficients for vertical scale factors at w-, t- levels 1187 !!gm bug : define it from analytical function, not like juste bellow.... 1188 !!gm or betteroffer the 2 possibilities.... 1189 DO jk = 1, jpkm1 1190 esigt(jk ) = gsigw(jk+1) - gsigw(jk) 1191 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1192 END DO 1193 esigw( 1 ) = esigw( 2 ) 1194 esigt(jpk) = esigt(jpkm1) 1195 1196 !!gm original form 1197 !!org DO jk = 1, jpk 1198 !!org esigt(jk)=fsdsig( FLOAT(jk) ) 1199 !!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 ) 1200 !!org END DO 1201 !!gm 1202 ! 1203 ! Coefficients for vertical depth as the sum of e3w scale factors 1506 1204 gsi3w(1) = 0.5 * esigw(1) 1507 1205 DO jk = 2, jpk 1508 gsi3w(jk) = gsi3w(jk-1)+ esigw(jk) 1509 END DO 1510 1511 !!bug gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1206 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) 1207 END DO 1208 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1512 1209 DO jk = 1, jpk 1513 !! change the solver stat.... 1514 !! zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1) 1515 !! gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft) 1516 gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*(FLOAT(jk)-0.5)/FLOAT(jpkm1)) 1517 !! zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 1518 !! gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw) 1519 !! gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw) 1520 gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1)) 1521 gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1)) 1522 END DO 1523 1524 !!bug gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1210 zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1) 1211 zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 1212 gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft) 1213 gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw) 1214 gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw) 1215 END DO 1216 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1525 1217 DO jj = 1, jpj 1526 1218 DO ji = 1, jpi … … 1530 1222 e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 1531 1223 e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1)) 1532 1224 ! 1533 1225 e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1)) 1534 1226 e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1)) … … 1537 1229 END DO 1538 1230 END DO 1539 1540 ! HYBRID 1231 ! 1232 ! HYBRID : 1541 1233 DO jj = 1, jpj 1542 1234 DO ji = 1, jpi … … 1547 1239 END DO 1548 1240 END DO 1549 1550 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1551 1552 1553 ! =========== 1554 ! Zoom domain 1555 ! =========== 1556 1557 IF( lzoom ) CALL zgr_bat_zoom 1558 1559 ! 2.4 Control print 1560 1561 IF(lwp) THEN 1241 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 1242 & ' MAX ', MAXVAL( mbathy(:,:) ) 1243 1244 1245 ! ! =========== 1246 IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain 1247 ! ! =========== 1248 1249 ! ! ============= 1250 IF(lwp) THEN ! Control print 1251 ! ! ============= 1562 1252 WRITE(numout,*) 1563 1253 WRITE(numout,*) ' domzgr: vertical coefficients for model level' 1564 WRITE(numout,9400) 1565 WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk) 1566 ENDIF 1567 9400 FORMAT(9x,' level gsigt gsigw esigt esigw gsi3w') 1568 9410 FORMAT(10x,i4,5f11.4) 1569 1570 IF(lwp) THEN 1254 WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" ) 1255 WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk ) 1256 ENDIF 1257 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 1258 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1259 WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), & 1260 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) ) 1261 WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), & 1262 & ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), & 1263 & ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), & 1264 & ' w ', MINVAL( fse3w (:,:,:) ) 1265 1266 WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), & 1267 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) ) 1268 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), & 1269 & ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), & 1270 & ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), & 1271 & ' w ', MAXVAL( fse3w (:,:,:) ) 1272 ENDIF 1273 ! 1274 IF(lwp) THEN ! selected vertical profiles 1571 1275 WRITE(numout,*) 1572 1276 WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 1573 1277 WRITE(numout,*) ' ~~~~~~ --------------------' 1574 WRITE(numout, 9420)1575 WRITE(numout, 9430) (jk,fsdept(1,1,jk),fsdepw(1,1,jk), &1576 fse3t (1,1,jk),fse3w (1,1,jk),jk=1,jpk)1278 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1279 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), & 1280 & fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk ) 1577 1281 DO jj = mj0(20), mj1(20) 1578 1282 DO ji = mi0(20), mi1(20) … … 1580 1284 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1581 1285 WRITE(numout,*) ' ~~~~~~ --------------------' 1582 WRITE(numout, 9420)1583 WRITE(numout, 9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk), &1584 & fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)1286 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1287 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1288 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1585 1289 END DO 1586 1290 END DO … … 1590 1294 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 1591 1295 WRITE(numout,*) ' ~~~~~~ --------------------' 1592 WRITE(numout,9420) 1593 WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk), & 1594 & fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk) 1595 END DO 1596 END DO 1597 ENDIF 1598 1599 9420 FORMAT(9x,' level gdept gdepw gde3w e3t e3w ') 1600 9430 FORMAT(10x,i4,4f9.2) 1601 1602 !!bug gm no more necessary? if ! defined key_helsinki 1296 WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')") 1297 WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), & 1298 & fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk ) 1299 END DO 1300 END DO 1301 ENDIF 1302 1303 !!gm bug? no more necessary? if ! defined key_helsinki 1603 1304 DO jk = 1, jpk 1604 1305 DO jj = 1, jpj 1605 1306 DO ji = 1, jpi 1606 1307 IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN 1607 IF(lwp) THEN 1608 WRITE(numout,*) 1609 WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 ' 1610 WRITE(numout,*) ' ========= --------------- ' 1611 WRITE(numout,*) 1612 WRITE(numout,*) ' point (i,j,k)= ',ji,jj,jk 1613 WRITE(numout,*) 1614 ENDIF 1615 STOP 'domzgr' 1308 WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1309 CALL ctl_stop( ctmp1 ) 1616 1310 ENDIF 1617 1311 IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN 1618 IF(lwp) THEN 1619 WRITE(numout,*) 1620 WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 ' 1621 WRITE(numout,*) ' ========= ------------------ ' 1622 WRITE(numout,*) 1623 WRITE(numout,*) ' point (i,j,k)= ', ji, jj, jk 1624 WRITE(numout,*) 1625 ENDIF 1626 STOP 'domzgr' 1312 WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1313 CALL ctl_stop( ctmp1 ) 1627 1314 ENDIF 1628 1315 END DO 1629 1316 END DO 1630 1317 END DO 1631 !! bug gm#endif1632 1318 !!gm bug #endif 1319 ! 1633 1320 END SUBROUTINE zgr_sco 1634 1321 1635 1322 #endif 1636 1637 1323 1638 1324 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
