Changeset 454 for trunk/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2006-05-10T18:47:31+02:00 (18 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90
r450 r454 11 11 !! zgr_bat_ctl : check the bathymetry files 12 12 !! zgr_z : reference z-coordinate 13 !! zgr_zco : z-coordinate 13 14 !! zgr_zps : z-coordinate with partial steps 14 !! zgr_s 15 !! zgr_sco : s-coordinate 15 16 !!--------------------------------------------------------------------- 16 17 !! * Modules used … … 30 31 PUBLIC dom_zgr ! called by dom_init.F90 31 32 33 !! * Module variables 34 REAL(wp) :: & !!: Namelist nam_zgr_sco 35 sbot_min = 300. , & !: minimum depth of s-bottom surface (>0) (m) 36 sbot_max = 5250. , & !: maximum depth of s-bottom surface (= ocean depth) (>0) (m) 37 theta = 6.0 , & !: surface control parameter (0<=theta<=20) 38 thetb = 0.75, & !: bottom control parameter (0<=thetb<= 1) 39 r_max = 0.15 !: maximum cut-off r-value allowed (0<r_max<1) 40 41 42 32 43 !! * Substitutions 33 44 # include "domzgr_substitute.h90" … … 35 46 !!---------------------------------------------------------------------- 36 47 !! OPA 9.0 , LOCEAN-IPSL (2005) 37 !! $Header$38 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt39 48 !!---------------------------------------------------------------------- 40 49 … … 49 58 !! 50 59 !! ** Method reference vertical coordinate 51 !! Z-coordinates : The depth of model levels is defined 52 !! from an analytical function the derivative of which gives 53 !! the vertical scale factors. 54 !! both depth and scale factors only depend on k (1d arrays). 55 !! w-level: gdepw = fsdep(k) 56 !! e3w(k) = dk(fsdep)(k) = fse3(k) 57 !! t-level: gdept = fsdep(k+0.5) 58 !! e3t(k) = dk(fsdep)(k+0.5) = fse3(k+0.5) 59 !! 60 !! ** Action : - gdept, gdepw : depth of T- and W-point (m) 61 !! - e3t, e3w : scale factors at T- and W-levels (m) 62 !! 63 !! Reference : 64 !! Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 60 !! 61 !! ** Action : 65 62 !! 66 63 !! History : 67 64 !! 9.0 ! 03-08 (G. Madec) original code 65 !! 9.0 ! 05-10 (A. Beckmann) modifications for hybrid s-ccordinates 68 66 !!---------------------------------------------------------------------- 69 67 INTEGER :: ioptio = 0 ! temporary integer 70 !!---------------------------------------------------------------------- 68 69 NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco 70 !!---------------------------------------------------------------------- 71 72 ! Read Namelist nam_zgr : vertical coordinate' 73 ! --------------------- 74 REWIND ( numnam ) 75 READ ( numnam, nam_zgr ) 76 77 ! Parameter control and print 78 ! --------------------------- 79 ! Control print 80 IF(lwp) THEN 81 WRITE(numout,*) 82 WRITE(numout,*) 'dom_zgr : vertical coordinate' 83 WRITE(numout,*) '~~~~~~~' 84 WRITE(numout,*) ' Namelist nam_zgr : set vertical coordinate' 85 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 86 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 87 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 88 ENDIF 71 89 72 90 ! Check Vertical coordinate options 73 ! ---------------------------------74 91 ioptio = 0 75 IF( lk_sco ) THEN 76 IF(lwp) WRITE(numout,*) 77 IF(lwp) WRITE(numout,*) 'dom_zgr : s-coordinate' 78 IF(lwp) WRITE(numout,*) '~~~~~~~' 79 ioptio = ioptio + 1 80 ENDIF 81 IF( lk_zps ) THEN 82 IF(lwp) WRITE(numout,*) 83 IF(lwp) WRITE(numout,*) 'dom_zgr : z-coordinate with partial steps' 84 IF(lwp) WRITE(numout,*) '~~~~~~~' 85 ioptio = ioptio + 1 86 ENDIF 87 IF( ioptio == 0 ) THEN 88 IF(lwp) WRITE(numout,*) 89 IF(lwp) WRITE(numout,*) 'dom_zgr : z-coordinate' 90 IF(lwp) WRITE(numout,*) '~~~~~~~' 91 ENDIF 92 93 IF ( ioptio > 1 ) THEN 92 IF( ln_zco ) ioptio = ioptio + 1 93 IF( ln_zps ) ioptio = ioptio + 1 94 IF( ln_sco ) ioptio = ioptio + 1 95 IF ( ioptio /= 1 ) THEN 94 96 IF(lwp) WRITE(numout,cform_err) 95 IF(lwp) WRITE(numout,*) ' several vertical coordinate options used'97 IF(lwp) WRITE(numout,*) ' none or several vertical coordinate options used' 96 98 nstop = nstop + 1 97 99 ENDIF 100 101 IF( ln_zco ) THEN 102 IF(lwp) WRITE(numout,*) ' z-coordinate with reduced incore memory requirement' 103 IF( ln_zps .OR. ln_sco ) THEN 104 IF(lwp) WRITE(numout,cform_err) 105 IF(lwp) WRITE(numout,*) ' reduced memory with zps or sco option is impossible' 106 nstop = nstop + 1 107 ENDIF 108 ENDIF 109 98 110 99 111 ! Build the vertical coordinate system 100 112 ! ------------------------------------ 101 113 102 CALL zgr_z ! Reference z-coordinate system 103 104 CALL zgr_bat ! Bathymetry fields (levels and meters) 105 106 CALL zgr_zps ! Partial step z-coordinate 107 108 CALL zgr_s ! s-coordinate 114 CALL zgr_z ! Reference z-coordinate system (always called) 115 116 CALL zgr_bat ! Bathymetry fields (levels and meters) 117 118 IF( ln_zco ) CALL zgr_zco ! z-coordinate 119 120 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 121 122 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 123 124 !!bug gm control print: 125 write(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 126 write(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), & 127 & ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) ) 128 write(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), & 129 & ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), & 130 & ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), & 131 & ' w ', MINVAL( fse3w(:,:,:) ) 132 133 write(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), & 134 & ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) ) 135 write(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), & 136 & ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), & 137 & ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), & 138 & ' w ', MAXVAL( fse3w(:,:,:) ) 139 !!!bug gm 109 140 110 141 END SUBROUTINE dom_zgr … … 122 153 !! function the derivative of which gives the scale factors. 123 154 !! both depth and scale factors only depend on k (1d arrays). 124 !! w-level: gdepw = fsdep(k)125 !! e3w (k) = dk(fsdep)(k) = fse3(k)126 !! t-level: gdept = fsdep(k+0.5)127 !! e3t (k) = dk(fsdep)(k+0.5) = fse3(k+0.5)128 !! 129 !! ** Action : - gdept , gdepw: depth of T- and W-point (m)130 !! - e3t , e3w: scale factors at T- and W-levels (m)155 !! w-level: gdepw_0 = fsdep(k) 156 !! e3w_0(k) = dk(fsdep)(k) = fse3(k) 157 !! t-level: gdept_0 = fsdep(k+0.5) 158 !! e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5) 159 !! 160 !! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m) 161 !! - e3t_0, e3w_0 : scale factors at T- and W-levels (m) 131 162 !! 132 163 !! Reference : … … 174 205 WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates' 175 206 WRITE(numout,*) ' ~~~~~~~' 176 IF 207 IF( ppkth == 0. ) THEN 177 208 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 178 209 WRITE(numout,*) ' Total depth :', zhmax 179 210 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 180 211 ELSE 181 IF 212 IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN 182 213 WRITE(numout,*) ' zsur, za0, za1 computed from ' 183 214 WRITE(numout,*) ' zdzmin = ', zdzmin … … 196 227 ! Reference z-coordinate (depth - scale factor at T- and W-points) 197 228 ! ====================== 198 IF ( ppkth == 0.) THEN ! uniform vertical grid199 200 za1 = zhmax /FLOAT(jpk-1)229 IF( ppkth == 0.e0 ) THEN ! uniform vertical grid 230 231 za1 = zhmax / FLOAT(jpk-1) 201 232 DO jk = 1, jpk 202 233 zw = FLOAT( jk ) 203 234 zt = FLOAT( jk ) + 0.5 204 gdepw (jk) = ( zw - 1 ) * za1205 gdept (jk) = ( zt - 1 ) * za1206 e3w (jk) = za1207 e3t (jk) = za1235 gdepw_0(jk) = ( zw - 1 ) * za1 236 gdept_0(jk) = ( zt - 1 ) * za1 237 e3w_0 (jk) = za1 238 e3t_0 (jk) = za1 208 239 END DO 209 240 … … 213 244 zw = FLOAT( jk ) 214 245 zt = FLOAT( jk ) + 0.5 215 gdepw (jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) ) )216 gdept (jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) ) )217 e3w (jk) = za0 + za1 * TANH( (zw-zkth)/zacr )218 e3t (jk) = za0 + za1 * TANH( (zt-zkth)/zacr )219 END DO 220 gdepw (1) = 0.e0! force first w-level to be exactly at zero246 gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) ) ) 247 gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) ) ) 248 e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth)/zacr ) 249 e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth)/zacr ) 250 END DO 251 gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero 221 252 222 253 ENDIF … … 229 260 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' 230 261 WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" ) 231 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept (jk), gdepw(jk), e3t(jk), e3w(jk), jk = 1, jpk )262 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk ) 232 263 ENDIF 233 264 234 265 DO jk = 1, jpk 235 IF( e3w (jk) <= 0. .OR. e3t(jk) <= 0. ) THEN266 IF( e3w_0(jk) <= 0. .OR. e3t_0(jk) <= 0. ) THEN 236 267 IF(lwp) WRITE(numout,cform_err) 237 268 IF(lwp) WRITE(numout,*) ' e3w or e3t =< 0 ' 238 269 nstop = nstop + 1 239 270 ENDIF 240 IF( gdepw (jk) < 0. .OR. gdept(jk) < 0.) THEN271 IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0.) THEN 241 272 IF(lwp) WRITE(numout,cform_err) 242 273 IF(lwp) WRITE(numout,*) ' gdepw or gdept < 0 ' … … 284 315 !! History : 285 316 !! 9.0 ! 03-08 (G. Madec) Original code 317 !! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates 286 318 !!---------------------------------------------------------------------- 287 319 !! * Modules used … … 295 327 INTEGER :: & 296 328 ipi, ipj, ipk, & ! temporary integers 297 itime, 329 itime, ih, & ! " " 298 330 ii_bump, ij_bump ! bump center position 299 331 INTEGER, DIMENSION (1) :: istep … … 302 334 REAL(wp) :: & 303 335 r_bump, h_bump, h_oce, & ! bump characteristics 304 zi, zj, zdate0, zdt 336 zi, zj, zdate0, zdt, zh ! temporary scalars 305 337 REAL(wp), DIMENSION(jpidta,jpjdta) :: & 306 zlamt, zphit, & ! temporaryworkspace (NetCDF read)338 zlamt, zphit, & ! 2D workspace (NetCDF read) 307 339 zdta ! global domain scalar data 308 340 REAL(wp), DIMENSION(jpk) :: & 309 zdept ! temporaryworkspace (NetCDF read)341 zdept ! 1D workspace (NetCDF read) 310 342 !!---------------------------------------------------------------------- 311 343 … … 327 359 328 360 idta(:,:) = jpkm1 ! flat basin 329 zdta(:,:) = gdepw (jpk)361 zdta(:,:) = gdepw_0(jpk) 330 362 331 363 ELSE ! bump … … 335 367 ii_bump = jpidta / 3 + 3 ! i-index of the bump center 336 368 ij_bump = jpjdta / 2 ! j-index of the bump center 337 r_bump = 6 ! bump radius (index)338 h_bump = 240.e0 ! bump height (meters)339 h_oce = gdepw (jpk)! background ocean depth (meters)369 r_bump = 6 ! bump radius (index) 370 h_bump = 240.e0 ! bump height (meters) 371 h_oce = gdepw_0(jpk) ! background ocean depth (meters) 340 372 IF(lwp) WRITE(numout,*) ' bump characteristics: ' 341 373 IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump … … 351 383 END DO 352 384 END DO 385 353 386 ! idta : 354 idta(:,:) = jpkm1 355 DO jk = 1, jpkm1 356 DO jj = 1, jpjdta 357 DO ji = 1, jpidta 358 IF( gdept(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept(jk+1) ) idta(ji,jj) = jk 387 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 388 idta(:,:) = jpkm1 389 ELSE ! z-coordinate (zco or zps): step-like topography 390 idta(:,:) = jpkm1 391 DO jk = 1, jpkm1 392 DO jj = 1, jpjdta 393 DO ji = 1, jpidta 394 IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk 395 END DO 359 396 END DO 360 397 END DO 361 END DO398 ENDIF 362 399 ENDIF 363 400 … … 372 409 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0.e0 373 410 ELSE 374 idta( : , 1 ) = 0 ; zdta( : , 1 ) = 0.e0 375 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0 376 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0.e0 377 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0.e0 411 ih = 0 ; zh = 0.e0 412 IF( ln_sco ) zh = jpkm1 ; IF( ln_sco ) zh = h_oce 413 idta( : , 1 ) = ih ; zdta( : , 1 ) = zh 414 idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh 415 idta( 1 , : ) = ih ; zdta( 1 , : ) = zh 416 idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh 378 417 ENDIF 379 418 380 419 ! EEL R5 configuration with east and west open boundaries. 381 420 ! Two rows of zeroes are needed at the south and north for OBCs 382 ! This is for compatibility with the rigid lid option.383 421 384 422 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN 385 idta( : , 2 ) = 0 ; zdta( : , 2 ) = 0.e0 386 idta( : ,jpjdta-1) = 0 ; zdta( : ,jpjdta-1) = 0.e0 423 ih = 0 ; zh = 0.e0 424 IF( ln_sco ) zh = jpkm1 ; IF( ln_sco ) zh = h_oce 425 idta( : , 2 ) = jpkm1 ; zdta( : , 2 ) = h_oce 426 idta( : ,jpjdta-1) = jpkm1 ; zdta( : ,jpjdta-1) = h_oce 387 427 ENDIF 388 428 … … 390 430 ELSEIF( ntopo == 1 ) THEN ! read in file ! 391 431 ! ! =============== ! 392 IF( lk_zco ) THEN 393 clname = 'bathy_level.nc' ! Level bathymetry 432 433 clname = 'bathy_level.nc' ! Level bathymetry 434 #if defined key_agrif 435 IF( .NOT. Agrif_Root() ) THEN 436 clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 437 ENDIF 438 #endif 439 INQUIRE( FILE=clname, EXIST=llbon ) 440 IF( llbon ) THEN 441 IF(lwp) WRITE(numout,*) 442 IF(lwp) WRITE(numout,*) ' read level bathymetry in ', clname 443 IF(lwp) WRITE(numout,*) 444 ipi = jpidta ; ipj = jpjdta 445 ipk = 1 ; itime = 1 ; zdt = rdt 446 CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE., & 447 & ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum ) 448 CALL flinget( inum, 'Bathy_level', jpidta, jpjdta, 1, & 449 & itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) ) 450 CALL flinclo( inum ) 451 idta(:,:) = zdta(:,:) 452 ELSE 453 IF( ln_zco ) THEN 454 IF(lwp) WRITE(numout,cform_err) 455 IF(lwp) WRITE(numout,*)' zgr_bat : unable to read the file ', clname 456 nstop = nstop + 1 457 ELSE 458 IF(lwp) WRITE(numout,*)' zgr_bat : bathy_level will be computed from bathy_meter' 459 idta(:,:) = jpkm1 ! initialisation 460 ENDIF 461 ENDIF 462 463 clname = 'bathy_meter.nc' ! meter bathymetry 394 464 #if defined key_agrif 395 465 IF( .NOT. Agrif_Root() ) THEN 396 466 clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 397 467 ENDIF 398 #endif 399 INQUIRE( FILE=clname, EXIST=llbon ) 400 IF( llbon ) THEN 401 IF(lwp) WRITE(numout,*) 402 IF(lwp) WRITE(numout,*) ' read level bathymetry in ', clname 403 IF(lwp) WRITE(numout,*) 404 itime = 1 405 ipi = jpidta 406 ipj = jpjdta 407 ipk = 1 408 zdt = rdt 409 CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE., & 410 ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum ) 411 CALL flinget( inum, 'Bathy_level', jpidta, jpjdta, 1, & 412 itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) ) 413 idta(:,:) = zdta(:,:) 414 CALL flinclo( inum ) 415 416 ELSE 417 IF(lwp) WRITE(numout,cform_err) 418 IF(lwp) WRITE(numout,*)' zgr_bat : unable to read the file', clname 419 nstop = nstop + 1 420 ENDIF 421 422 ELSEIF( lk_zps ) THEN 423 clname = 'bathy_meter.nc' ! meter bathymetry 424 #if defined key_agrif 425 IF( .NOT. Agrif_Root() ) THEN 426 clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname) 427 ENDIF 428 #endif 429 INQUIRE( FILE=clname, EXIST=llbon ) 430 IF( llbon ) THEN 431 IF(lwp) WRITE(numout,*) 432 IF(lwp) WRITE(numout,*) ' read meter bathymetry in ', clname 433 IF(lwp) WRITE(numout,*) 434 itime = 1 435 ipi = jpidta 436 ipj = jpjdta 437 ipk = 1 438 zdt = rdt 439 CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE., & 440 ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum ) 441 CALL flinget( inum, 'Bathymetry', jpidta, jpjdta, 1, & 442 itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) ) 443 CALL flinclo( inum ) 444 idta(:,:) = jpkm1 ! initialisation 445 ELSE 468 #endif 469 INQUIRE( FILE=clname, EXIST=llbon ) 470 IF( llbon ) THEN 471 IF(lwp) WRITE(numout,*) 472 IF(lwp) WRITE(numout,*) ' read meter bathymetry in ', clname 473 IF(lwp) WRITE(numout,*) 474 ipi = jpidta ; ipj = jpjdta 475 ipk = 1 ; itime = 1 ; zdt = rdt 476 CALL flinopen( clname, 1, jpidta, 1, jpjdta, .FALSE., & 477 & ipi, ipj, ipk, zlamt, zphit, zdept, itime, istep, zdate0, zdt, inum ) 478 CALL flinget( inum, 'Bathymetry', jpidta, jpjdta, 1, & 479 & itime, 1, 1, 1, jpidta, 1, jpjdta, zdta(:,:) ) 480 CALL flinclo( inum ) 481 ELSE 482 IF( ln_zps .OR. ln_sco ) THEN 446 483 IF(lwp) WRITE(numout,cform_err) 447 484 IF(lwp) WRITE(numout,*)' zgr_bat : unable to read the file', clname 448 485 nstop = nstop + 1 486 ELSE 487 zdta(:,:) = 0.e0 488 IF(lwp) WRITE(numout,*)' zgr_bat : bathy_meter not found, but not used, bathy array set to zero' 449 489 ENDIF 450 490 ENDIF … … 472 512 END DO 473 513 514 write(numout,*) ' MIN val mbathy 2 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 515 516 474 517 ! ======================= 475 518 ! NO closed seas or lakes … … 482 525 mbathy(ji,jj) = 0 ! suppress closed seas 483 526 bathy (ji,jj) = 0.e0 ! and lakes 527 !!bug gm bathy (ji,jj) = 300.e0 ! and lakes 484 528 END DO 485 529 END DO … … 585 629 !! History : 586 630 !! 9.0 ! 03-08 (G. Madec) Original code 631 !! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates 587 632 !!---------------------------------------------------------------------- 588 633 !! * Local declarations … … 621 666 DO ji = 2, jpim1 622 667 ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), & 623 mbathy(ji,jj-1),mbathy(ji,jj+1) )668 & mbathy(ji,jj-1),mbathy(ji,jj+1) ) 624 669 IF( ibtest < mbathy(ji,jj) ) THEN 625 670 IF(lwp) WRITE(numout,*) ' the number of ocean level at ', & 626 'grid-point (i,j) = ',ji,jj,' is changed from ', & 627 mbathy(ji,jj),' to ', ibtest 671 & 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest 628 672 mbathy(ji,jj) = ibtest 629 673 icompt = icompt + 1 … … 657 701 ENDIF 658 702 ELSE 659 mbathy( 1 ,:) = 0 660 mbathy(jpi,:) = 0 703 IF( ln_zco .OR. ln_zps ) THEN 704 mbathy( 1 ,:) = 0 705 mbathy(jpi,:) = 0 706 ELSE 707 mbathy( 1 ,:) = jpkm1 708 mbathy(jpi,:) = jpkm1 709 ENDIF 661 710 ENDIF 662 711 ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN … … 684 733 ! Boundary condition on mbathy 685 734 IF( .NOT.lk_mpp ) THEN 686 687 735 !!bug ??? y reflechir! 688 736 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab 689 690 737 zbathy(:,:) = FLOAT( mbathy(:,:) ) 691 738 CALL lbc_lnk( zbathy, 'T', 1. ) … … 729 776 730 777 778 SUBROUTINE zgr_zco 779 !!---------------------------------------------------------------------- 780 !! *** ROUTINE zgr_zco *** 781 !! 782 !! ** Purpose : define the z-coordinate system 783 !! 784 !! ** Method : set 3D coord. arrays to reference 1D array if lk_zco=T 785 !! 786 !! History : 787 !! ! 06-04 (R. Benshila, G. Madec) Original code 788 !!---------------------------------------------------------------------- 789 INTEGER :: jk 790 !!---------------------------------------------------------------------- 791 792 IF( .NOT.lk_zco ) THEN 793 DO jk = 1, jpk 794 fsdept(:,:,jk) = gdept_0(jk) 795 fsdepw(:,:,jk) = gdepw_0(jk) 796 fsde3w(:,:,jk) = gdepw_0(jk) 797 fse3t (:,:,jk) = e3t_0(jk) 798 fse3u (:,:,jk) = e3t_0(jk) 799 fse3v (:,:,jk) = e3t_0(jk) 800 fse3f (:,:,jk) = e3t_0(jk) 801 fse3w (:,:,jk) = e3w_0(jk) 802 fse3uw(:,:,jk) = e3w_0(jk) 803 fse3vw(:,:,jk) = e3w_0(jk) 804 END DO 805 ENDIF 806 807 END SUBROUTINE zgr_zco 808 809 810 731 811 # include "domzgr_zps.h90" 732 812 733 813 734 # include "domzgr_s.h90" 814 #if ! defined key_zco 815 !!---------------------------------------------------------------------- 816 !! NOT 'key_zco' : 3D gdep. and e3. 817 !!---------------------------------------------------------------------- 818 819 SUBROUTINE zgr_sco 820 !!---------------------------------------------------------------------- 821 !! *** ROUTINE zgr_sco *** 822 !! 823 !! ** Purpose : define the s-coordinate system 824 !! 825 !! ** Method : s-coordinate 826 !! The depth of model levels is defined as the product of an 827 !! analytical function by the local bathymetry, while the vertical 828 !! scale factors are defined as the product of the first derivative 829 !! of the analytical function by the bathymetry. 830 !! (this solution save memory as depth and scale factors are not 831 !! 3d fields) 832 !! - Read bathymetry (in meters) at t-point and compute the 833 !! bathymetry at u-, v-, and f-points. 834 !! hbatu = mi( hbatt ) 835 !! hbatv = mj( hbatt ) 836 !! hbatf = mi( mj( hbatt ) ) 837 !! - Compute gsigt, gsigw, esigt, esigw from an analytical 838 !! function and its derivative given as statement function. 839 !! gsigt(k) = fssig (k+0.5) 840 !! gsigw(k) = fssig (k ) 841 !! esigt(k) = fsdsig(k+0.5) 842 !! esigw(k) = fsdsig(k ) 843 !! This routine is given as an example, it must be modified 844 !! following the user s desiderata. nevertheless, the output as 845 !! well as the way to compute the model levels and scale factors 846 !! must be respected in order to insure second order a!!uracy 847 !! schemes. 848 !! 849 !! Reference : 850 !! Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 851 !! 852 !! History : 853 !! ! 95-12 (G. Madec) Original code : s vertical coordinate 854 !! ! 97-07 (G. Madec) lbc_lnk call 855 !! ! 97-04 (J.-O. Beismann) 856 !! 8.5 ! 02-06 (G. Madec) F90: Free form and module 857 !! 9.0 ! 05-10 (A. Beckmann) new stretching function 858 !!---------------------------------------------------------------------- 859 !! * Local declarations 860 INTEGER :: ji, jj, jk, jl 861 REAL(wp) :: fssig, fsdsig, pfkk 862 863 INTEGER :: iip1, ijp1, iim1, ijm1 864 REAL(wp) :: & 865 fssigt, fssigw, zcoeft, zcoefw, & 866 zrmax, ztaper 867 868 REAL(wp), DIMENSION(jpi,jpj) :: & 869 zenv, ztmp, zmsk, zri, zrj, zhbat 870 871 NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max 872 !!---------------------------------------------------------------------- 873 874 ! a. Sigma stretching coefficient 875 fssigt(pfkk) = (tanh(theta*(-(pfkk-0.5)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) & 876 * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta)) 877 fssigw(pfkk) = (tanh(theta*(-(pfkk-1.0)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) & 878 * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta)) 879 880 ! b. Vertical derivative of sigma stretching coefficient 881 882 !!bug fsdsig(pfkk)= put here the analytical derivative of fssigt and w ... 883 !!---------------------------------------------------------------------- 884 885 ! Read Namelist nam_zgr_sco : sigma-stretching parameters 886 ! ------------------------- 887 REWIND( numnam ) 888 READ ( numnam, nam_zgr_sco ) 889 890 IF(lwp) THEN 891 WRITE(numout,*) 892 WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate' 893 WRITE(numout,*) '~~~~~~~~~~~' 894 WRITE(numout,*) ' Namelist nam_zgr_sco' 895 WRITE(numout,*) ' sigma-stretching coeffs ' 896 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max 897 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min 898 WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta 899 WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb 900 WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max 901 ENDIF 902 903 ! ??? 904 ! ------------------ 905 hift(:,:) = sbot_min 906 hifu(:,:) = sbot_min 907 hifv(:,:) = sbot_min 908 hiff(:,:) = sbot_min 909 910 911 ! set maximum ocean depth 912 ! ----------------------- 913 DO jj = 1, jpj 914 DO ji = 1, jpi 915 bathy(ji,jj) = MIN( sbot_max, bathy(ji,jj) ) 916 END DO 917 END DO 918 919 ! Define the envelop bathymetry 920 ! ============================= 921 ! Smooth the bathymetry (if required) 922 923 scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea) 924 scobot(:,:) = bathy(:,:) ! ocean bottom depth 925 926 927 ! use r-value to create hybrid coordinates 928 929 DO jj = 1, jpj 930 DO ji = 1, jpi 931 zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min ) 932 END DO 933 END DO 934 935 jl = 0 936 zrmax = 1.e0 937 ! ! ================ ! 938 DO WHILE ( jl <= 10000 .AND. zrmax > r_max ) ! Iterative loop ! 939 ! ! ================ ! 940 jl = jl + 1 941 942 zrmax = 0.e0 943 zmsk(:,:) = 0.e0 944 945 DO jj = 1, nlcj 946 DO ji = 1, nlci 947 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 948 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 949 zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 950 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 951 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 952 IF( zri(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0 953 IF( zri(ji,jj) > r_max ) zmsk(iip1,jj ) = 1.0 954 IF( zrj(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0 955 IF( zrj(ji,jj) > r_max ) zmsk(ji ,ijp1) = 1.0 956 END DO 957 END DO 958 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 959 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 960 ztmp(:,:) = zmsk(:,:) 961 CALL lbc_lnk( zmsk, 'T', 1. ) 962 DO jj = 1, nlcj 963 DO ji = 1, nlci 964 zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) ) 965 END DO 966 END DO 967 968 WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) ) 969 970 DO jj = 1, nlcj 971 DO ji = 1, nlci 972 iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci) 973 ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj) 974 iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) 975 ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) 976 IF( zmsk(ji,jj) == 1.0 ) THEN 977 ztmp(ji,jj) = ( & 978 & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & 979 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) & 980 & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & 981 & ) / ( & 982 & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & 983 & + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) & 984 & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & 985 & ) 986 ENDIF 987 END DO 988 END DO 989 990 DO jj = 1, nlcj 991 DO ji = 1, nlci 992 IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 993 END DO 994 END DO 995 996 ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero 997 ztmp(:,:) = zenv(:,:) 998 CALL lbc_lnk( zenv, 'T', 1. ) 999 DO jj = 1, nlcj 1000 DO ji = 1, nlci 1001 IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj) 1002 END DO 1003 END DO 1004 ! ! ================ ! 1005 END DO ! End loop ! 1006 ! ! ================ ! 1007 1008 ! save envelop bathymetry in hbatt 1009 hbatt(:,:) = zenv(:,:) 1010 1011 !!bug gm : CAUTION: tapering hard coded !!!! orca2 only 1012 !!bug gm NOT valid in mpp ===> taper have been changed 1013 1014 DO jj = 1, jpj 1015 DO ji = 1, jpi 1016 !!bug mpp ztaper = EXP( -(FLOAT(jj-74)/10.)**2 ) 1017 ztaper = EXP( -(gphit(ji,jj)/8)**2 ) 1018 hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper) 1019 END DO 1020 END DO 1021 1022 DO jj = 1, jpj 1023 zrmax = EXP( -(gphit(10,jj)/8)**2 ) 1024 ztaper = EXP( -(FLOAT(jj-74)/10.)**2 ) 1025 write(numout,*) jj, (FLOAT(jj-74)/10.), ztaper,(gphit(10,jj)/8), zrmax 1026 END DO 1027 1028 write(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) 1029 write(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) ) 1030 1031 ! Control print 1032 IF(lwp) THEN 1033 WRITE(numout,*) 1034 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 1035 WRITE(numout,*) 1036 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout ) 1037 ENDIF 1038 1039 ! 1. computation of hbatu, hbatv, hbatf fields 1040 ! -------------------------------------------- 1041 1042 hbatu(:,:) = sbot_min 1043 hbatv(:,:) = sbot_min 1044 hbatf(:,:) = sbot_min 1045 1046 IF(lwp) THEN 1047 WRITE(numout,*) 1048 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min 1049 WRITE(numout,*) 1050 ENDIF 1051 1052 DO jj = 1, jpjm1 1053 DO ji = 1, jpim1 1054 hbatu(ji,jj)= 0.5 *( hbatt(ji ,jj)+hbatt(ji+1,jj ) ) 1055 hbatv(ji,jj)= 0.5 *( hbatt(ji ,jj)+hbatt(ji ,jj+1) ) 1056 hbatf(ji,jj)= 0.25*( hbatt(ji ,jj)+hbatt(ji ,jj+1) & 1057 +hbatt(ji+1,jj)+hbatt(ji+1,jj+1) ) 1058 END DO 1059 END DO 1060 1061 ! Apply lateral boundary condition 1062 ! CAUTION: retain non zero value in the initial file 1063 ! this should be OK for orca cfg, not for EEL 1064 zhbat(:,:) = hbatu(:,:) 1065 CALL lbc_lnk( hbatu, 'U', 1. ) 1066 DO jj = 1, jpj 1067 DO ji = 1, jpi 1068 IF( hbatu(ji,jj) == 0.e0 ) THEN 1069 IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = sbot_min 1070 IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj) 1071 ENDIF 1072 END DO 1073 END DO 1074 1075 zhbat(:,:) = hbatv(:,:) 1076 CALL lbc_lnk( hbatv, 'V', 1. ) 1077 DO jj = 1, jpj 1078 DO ji = 1, jpi 1079 IF( hbatv(ji,jj) == 0.e0 ) THEN 1080 IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = sbot_min 1081 IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj) 1082 ENDIF 1083 END DO 1084 END DO 1085 1086 zhbat(:,:) = hbatf(:,:) 1087 CALL lbc_lnk( hbatf, 'F', 1. ) 1088 DO jj = 1, jpj 1089 DO ji = 1, jpi 1090 IF( hbatf(ji,jj) == 0.e0 ) THEN 1091 IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = sbot_min 1092 IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj) 1093 ENDIF 1094 END DO 1095 END DO 1096 1097 1098 !!bug: key_helsinki a verifer 1099 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 1100 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 1101 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 1102 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 1103 1104 write(numout,*) ' MAX val hif t ', MAXVAL( hift(:,:) ), ' f ', MAXVAL( hiff(:,:) ), & 1105 & ' u ', MAXVAL( hifu(:,:) ), ' v ', MAXVAL( hifv(:,:) ) 1106 write(numout,*) ' MIN val hif t ', MINVAL( hift(:,:) ), ' f ', MINVAL( hiff(:,:) ), & 1107 & ' u ', MINVAL( hifu(:,:) ), ' v ', MINVAL( hifv(:,:) ) 1108 write(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & 1109 & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 1110 write(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & 1111 & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 1112 !! helsinki 1113 1114 ! 2. Computation of gsig and esig fields 1115 ! -------------------------------------- 1116 1117 ! 2.1 Coefficients for model level depth at w- and t-levels 1118 1119 !!bug gm : change the def of statment function.... 1120 DO jk = 1, jpk 1121 gsigw(jk) = -fssigw(FLOAT(jk)) 1122 gsigt(jk) = -fssigt(FLOAT(jk)) 1123 END DO 1124 1125 write(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1126 1127 !!org DO jk = 1, jpk 1128 !!org gsigw(jk) = -fssig ( FLOAT(jk) ) 1129 !!org gsigt(jk) = -fssig ( FLOAT(jk)+0.5) 1130 !!org END DO 1131 1132 ! 2.2 Coefficients for vertical scale factors at w-, t- levels 1133 1134 !!bug gm : define it from analytical function, not like juste bellow.... 1135 !! or betteroffer the 2 possibilities.... 1136 DO jk=2,jpk 1137 esigw(jk)=gsigt(jk)-gsigt(jk-1) 1138 END DO 1139 DO jk=1,jpkm1 1140 esigt(jk)=gsigw(jk+1)-gsigw(jk) 1141 END DO 1142 esigw(1)=esigw(2) 1143 esigt(jpk)=esigt(jpkm1) 1144 1145 !!org DO jk = 1, jpk 1146 !!org esigw(jk)=fsdsig( FLOAT(jk)) 1147 !!org esigt(jk)=fsdsig( FLOAT(jk)+0.5) 1148 !!org END DO 1149 1150 ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors 1151 1152 gsi3w(1) = 0.5 * esigw(1) 1153 DO jk = 2, jpk 1154 gsi3w(jk) = gsi3w(jk-1)+ esigw(jk) 1155 END DO 1156 1157 !!bug gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1158 DO jk = 1, jpk 1159 !! change the solver stat.... 1160 !! zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1) 1161 !! gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft) 1162 gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*(FLOAT(jk)-0.5)/FLOAT(jpkm1)) 1163 !! zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1) 1164 !! gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw) 1165 !! gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoefw) 1166 gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1)) 1167 gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*FLOAT(jk-1)/FLOAT(jpkm1)) 1168 END DO 1169 1170 !!bug gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 1171 DO jj = 1, jpj 1172 DO ji = 1, jpi 1173 DO jk = 1, jpk 1174 e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1)) 1175 e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 1176 e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 1177 e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1)) 1178 1179 e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1)) 1180 e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1)) 1181 e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1)) 1182 END DO 1183 END DO 1184 END DO 1185 1186 ! HYBRID 1187 DO jj = 1, jpj 1188 DO ji = 1, jpi 1189 DO jk = 1, jpkm1 1190 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1191 IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0 1192 END DO 1193 END DO 1194 END DO 1195 1196 write(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 1197 1198 1199 ! =========== 1200 ! Zoom domain 1201 ! =========== 1202 1203 IF( lzoom ) CALL zgr_bat_zoom 1204 1205 ! 2.4 Control print 1206 1207 IF(lwp) THEN 1208 WRITE(numout,*) 1209 WRITE(numout,*) ' domzgr: vertical coefficients for model level' 1210 WRITE(numout,9400) 1211 WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk) 1212 ENDIF 1213 9400 FORMAT(9x,' level gsigt gsigw esigt esigw gsi3w') 1214 9410 FORMAT(10x,i4,5f11.4) 1215 1216 IF(lwp) THEN 1217 WRITE(numout,*) 1218 WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 1219 WRITE(numout,*) ' ~~~~~~ --------------------' 1220 WRITE(numout,9420) 1221 WRITE(numout,9430) (jk,fsdept(1,1,jk),fsdepw(1,1,jk), & 1222 fse3t (1,1,jk),fse3w (1,1,jk),jk=1,jpk) 1223 WRITE(numout,*) 1224 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(20,20), hbatt(20,20) 1225 WRITE(numout,*) ' ~~~~~~ --------------------' 1226 WRITE(numout,9420) 1227 WRITE(numout,9430) (jk,fsdept(20,20,jk),fsdepw(20,20,jk), & 1228 fse3t (20,20,jk),fse3w (20,20,jk),jk=1,jpk) 1229 WRITE(numout,*) 1230 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(100,74), hbatt(100,74) 1231 WRITE(numout,*) ' ~~~~~~ --------------------' 1232 WRITE(numout,9420) 1233 WRITE(numout,9430) (jk,fsdept(100,74,jk),fsdepw(100,74,jk), & 1234 fse3t (100,74,jk),fse3w (100,74,jk),jk=1,jpk) 1235 ENDIF 1236 1237 9420 FORMAT(9x,' level gdept gdepw gde3w e3t e3w ') 1238 9430 FORMAT(10x,i4,4f9.2) 1239 1240 !!bug gm no more necessary? if ! defined key_helsinki 1241 DO jk = 1, jpk 1242 DO jj = 1, jpj 1243 DO ji = 1, jpi 1244 IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN 1245 IF(lwp) THEN 1246 WRITE(numout,*) 1247 WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 ' 1248 WRITE(numout,*) ' ========= --------------- ' 1249 WRITE(numout,*) 1250 WRITE(numout,*) ' point (i,j,k)= ',ji,jj,jk 1251 WRITE(numout,*) 1252 ENDIF 1253 STOP 'domzgr' 1254 ENDIF 1255 IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN 1256 IF(lwp) THEN 1257 WRITE(numout,*) 1258 WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 ' 1259 WRITE(numout,*) ' ========= ------------------ ' 1260 WRITE(numout,*) 1261 WRITE(numout,*) ' point (i,j,k)= ', ji, jj, jk 1262 WRITE(numout,*) 1263 ENDIF 1264 STOP 'domzgr' 1265 ENDIF 1266 END DO 1267 END DO 1268 END DO 1269 !!bug gm #endif 1270 1271 END SUBROUTINE zgr_sco 1272 1273 #else 1274 !!---------------------------------------------------------------------- 1275 !! Default option : Empty routine 1276 !!---------------------------------------------------------------------- 1277 SUBROUTINE zgr_sco 1278 END SUBROUTINE zgr_sco 1279 #endif 735 1280 736 1281
Note: See TracChangeset
for help on using the changeset viewer.
