Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Property svn:eol-style deleted
r2465 r2528 7 7 !! ! 1997-07 (G. Madec) lbc_lnk call 8 8 !! ! 1997-04 (J.-O. Beismann) 9 !! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module10 !! - ! 2002-09 (A. de Miranda) rigid-lid + islands9 !! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module 10 !! - ! 2002-09 (A. de Miranda) rigid-lid + islands 11 11 !! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module 12 12 !! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function … … 14 14 !! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style 15 15 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 16 !! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level 16 17 !!---------------------------------------------------------------------- 17 18 … … 21 22 !! zgr_bat_zoom : modify the bathymetry field if zoom domain 22 23 !! zgr_bat_ctl : check the bathymetry files 24 !! zgr_bot_level: deepest ocean level for t-, u, and v-points 23 25 !! zgr_z : reference z-coordinate 24 26 !! zgr_zco : z-coordinate … … 28 30 !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!) 29 31 !!--------------------------------------------------------------------- 30 USE oce ! ocean dynamics and tracers 31 USE dom_oce ! ocean space and time domain 32 USE in_out_manager ! I/O manager 33 USE lib_mpp ! distributed memory computing library 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 USE closea ! closed seas 36 USE c1d ! 1D configuration 32 USE oce ! ocean variables 33 USE dom_oce ! ocean domain 34 USE closea ! closed seas 35 USE c1d ! 1D vertical configuration 36 USE in_out_manager ! I/O manager 37 USE iom ! I/O library 38 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 39 USE lib_mpp ! distributed memory computing library 37 40 38 41 IMPLICIT NONE 39 42 PRIVATE 40 43 41 PUBLIC dom_zgr ! called by dom_init.F90 42 43 !!gm DOCTOR name for the namelist parameter... 44 ! !!! ** Namelist namzgr_sco ** 45 REAL(wp) :: rn_sbot_min = 300. ! minimum depth of s-bottom surface (>0) (m) 46 REAL(wp) :: rn_sbot_max = 5250. ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 47 REAL(wp) :: rn_theta = 6.0 ! surface control parameter (0<=rn_theta<=20) 48 REAL(wp) :: rn_thetb = 0.75 ! bottom control parameter (0<=rn_thetb<= 1) 49 REAL(wp) :: rn_rmax = 0.15 ! maximum cut-off r-value allowed (0<rn_rmax<1) 50 LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 51 REAL(wp) :: rn_bb = 0.8 ! stretching parameter for song and haidvogel stretching 52 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 53 REAL(wp) :: rn_hc = 150. ! Critical depth for s-sigma coordinates 44 PUBLIC dom_zgr ! called by dom_init.F90 45 46 ! !!* Namelist namzgr_sco * 47 REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m) 48 REAL(wp) :: rn_sbot_max = 5250._wp ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 49 REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20) 50 REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1) 51 REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1) 52 LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T) 53 REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter for song and haidvogel stretching 54 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom) 55 REAL(wp) :: rn_hc = 150._wp ! Critical depth for s-sigma coordinates 54 56 55 57 !! * Substitutions … … 57 59 # include "vectopt_loop_substitute.h90" 58 60 !!---------------------------------------------------------------------- 59 !! NEMO/OPA 3. 0 , LOCEAN-IPSL (2008)61 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 60 62 !! $Id$ 61 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)63 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 62 64 !!---------------------------------------------------------------------- 63 64 65 CONTAINS 65 66 … … 75 76 !! - vertical coordinate (gdep., e3.) depending on the 76 77 !! coordinate chosen : 77 !! ln_zco=T z-coordinate (forced if lk_zco)78 !! ln_zco=T z-coordinate 78 79 !! ln_zps=T z-coordinate with partial steps 79 80 !! ln_zco=T s-coordinate … … 82 83 !!---------------------------------------------------------------------- 83 84 INTEGER :: ioptio = 0 ! temporary integer 84 ! !85 ! 85 86 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 86 87 !!---------------------------------------------------------------------- 87 88 88 REWIND ( numnam )! Read Namelist namzgr : vertical coordinate'89 READ 89 REWIND( numnam ) ! Read Namelist namzgr : vertical coordinate' 90 READ ( numnam, namzgr ) 90 91 91 92 IF(lwp) THEN ! Control print … … 103 104 IF( ln_zps ) ioptio = ioptio + 1 104 105 IF( ln_sco ) ioptio = ioptio + 1 105 IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 106 IF( lk_zco ) THEN 107 IF(lwp) WRITE(numout,*) ' z-coordinate with reduced incore memory requirement' 108 IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' ) 109 ENDIF 110 106 IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 107 ! 111 108 ! Build the vertical coordinate system 112 109 ! ------------------------------------ 113 CALL zgr_z! Reference z-coordinate system (always called)114 CALL zgr_bat! Bathymetry fields (levels and meters)115 IF( ln_zco ) CALL zgr_zco! z-coordinate116 IF( ln_zps ) CALL zgr_zps! Partial step z-coordinate117 IF( ln_sco ) CALL zgr_sco! s-coordinate or hybrid z-s coordinate118 ! 119 ! final adjustment of mbathy & check120 ! -----------------------------------121 IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain122 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isoated ocean points123 124 125 !!bug gm control print: 110 CALL zgr_z ! Reference z-coordinate system (always called) 111 CALL zgr_bat ! Bathymetry fields (levels and meters) 112 IF( ln_zco ) CALL zgr_zco ! z-coordinate 113 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 114 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 115 ! 116 ! final adjustment of mbathy & check 117 ! ----------------------------------- 118 IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain 119 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isoated ocean points 120 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 121 ! 122 ! 126 123 IF( nprint == 1 .AND. lwp ) THEN 127 124 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) … … 140 137 & ' w ', MAXVAL( fse3w(:,:,:) ) 141 138 ENDIF 142 !!!bug gm 143 139 ! 144 140 END SUBROUTINE dom_zgr 145 141 … … 171 167 REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90 172 168 REAL(wp) :: zrefdep ! depth of the reference level (~10m) 169 REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters 173 170 !!---------------------------------------------------------------------- 174 171 … … 177 174 zkth = ppkth ; zacr = ppacr 178 175 zdzmin = ppdzmin ; zhmax = pphmax 176 zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters 179 177 180 178 ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed 181 179 ! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 182 !183 180 IF( ppa1 == pp_to_be_computed .AND. & 184 181 & ppa0 == pp_to_be_computed .AND. & … … 192 189 ELSE 193 190 za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur 191 za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter 194 192 ENDIF 195 193 … … 198 196 WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates' 199 197 WRITE(numout,*) ' ~~~~~~~' 200 IF( ppkth == 0. ) THEN198 IF( ppkth == 0._wp ) THEN 201 199 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 202 200 WRITE(numout,*) ' Total depth :', zhmax 203 201 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 204 202 ELSE 205 IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0.) THEN203 IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN 206 204 WRITE(numout,*) ' zsur, za0, za1 computed from ' 207 205 WRITE(numout,*) ' zdzmin = ', zdzmin … … 214 212 WRITE(numout,*) ' zkth = ', zkth 215 213 WRITE(numout,*) ' zacr = ', zacr 214 IF( ldbletanh ) THEN 215 WRITE(numout,*) ' (Double tanh za2 = ', za2 216 WRITE(numout,*) ' parameters) zkth2= ', zkth2 217 WRITE(numout,*) ' zacr2= ', zacr2 218 ENDIF 216 219 ENDIF 217 220 ENDIF … … 220 223 ! Reference z-coordinate (depth - scale factor at T- and W-points) 221 224 ! ====================== 222 IF( ppkth == 0. e0) THEN ! uniform vertical grid225 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 223 226 za1 = zhmax / FLOAT(jpk-1) 224 227 DO jk = 1, jpk 225 228 zw = FLOAT( jk ) 226 zt = FLOAT( jk ) + 0.5 229 zt = FLOAT( jk ) + 0.5_wp 227 230 gdepw_0(jk) = ( zw - 1 ) * za1 228 231 gdept_0(jk) = ( zt - 1 ) * za1 … … 231 234 END DO 232 235 ELSE ! Madec & Imbard 1996 function 233 DO jk = 1, jpk 234 zw = FLOAT( jk ) 235 zt = FLOAT( jk ) + 0.5 236 gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) ) 237 gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) ) 238 e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth) / zacr ) 239 e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth) / zacr ) 240 END DO 241 gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero 236 IF( .NOT. ldbletanh ) THEN 237 DO jk = 1, jpk 238 zw = REAL( jk , wp ) 239 zt = REAL( jk , wp ) + 0.5_wp 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 ELSE 246 DO jk = 1, jpk 247 zw = FLOAT( jk ) 248 zt = FLOAT( jk ) + 0.5_wp 249 ! Double tanh function 250 gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) & 251 & + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) ) 252 gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) & 253 & + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) ) 254 e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) & 255 & + za2 * TANH( (zw-zkth2) / zacr2 ) 256 e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) & 257 & + za2 * TANH( (zt-zkth2) / zacr2 ) 258 END DO 259 ENDIF 260 gdepw_0(1) = 0._wp ! force first w-level to be exactly at zero 242 261 ENDIF 243 262 244 263 !!gm BUG in s-coordinate this does not work! 245 ! deepest/shallowest W level Above/Bel low ~10m246 zrefdep = 10. - ( 0.1*MINVAL(e3w_0) )! ref. depth with tolerance (10% of minimum layer thickness)247 nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Bel low ~10m248 nla10 = nlb10 - 1 ! deepest W level Above 264 ! deepest/shallowest W level Above/Below ~10m 265 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_0 ) ! ref. depth with tolerance (10% of minimum layer thickness) 266 nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 267 nla10 = nlb10 - 1 ! deepest W level Above ~10m 249 268 !!gm end bug 250 269 … … 256 275 ENDIF 257 276 DO jk = 1, jpk ! control positivity 258 IF( e3w_0 (jk) <= 0. e0 .OR. e3t_0 (jk) <= 0.e0 ) CALL ctl_stop( 'e3w or e3t =< 0 ' )259 IF( gdepw_0(jk) < 0. e0 .OR. gdept_0(jk) < 0.e0 ) CALL ctl_stop( 'gdepw or gdept < 0 ' )277 IF( e3w_0 (jk) <= 0._wp .OR. e3t_0 (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 ' ) 278 IF( gdepw_0(jk) < 0._wp .OR. gdept_0(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' ) 260 279 END DO 261 280 ! … … 289 308 !! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file 290 309 !! bathy is read in 'bathy_meter.nc' NetCDF file 291 !! C A U T I O N : mbathy will be modified during the initializa-292 !! tion phase to become the number of non-zero w-levels of a water293 !! column, with a minimum value of 1.294 310 !! 295 311 !! ** Action : - mbathy: level bathymetry (in level index) 296 312 !! - bathy : meter bathymetry (in meters) 297 313 !!---------------------------------------------------------------------- 298 USE iom299 !!300 314 INTEGER :: ji, jj, jl, jk ! dummy loop indices 301 315 INTEGER :: inum ! temporary logical unit 302 316 INTEGER :: ii_bump, ij_bump, ih ! bump center position 303 INTEGER :: ii0, ii1, ij0, ij1 !indices317 INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices 304 318 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 305 REAL(wp) :: zi , zj , zh ! temporaryscalars319 REAL(wp) :: zi, zj, zh, zhmin ! local scalars 306 320 INTEGER , DIMENSION(jpidta,jpjdta) :: idta ! global domain integer data 307 321 REAL(wp), DIMENSION(jpidta,jpjdta) :: zdta ! global domain scalar data … … 328 342 ii_bump = jpidta / 2 ! i-index of the bump center 329 343 ij_bump = jpjdta / 2 ! j-index of the bump center 330 r_bump = 50000. e0! bump radius (meters)331 h_bump = 2700.e0! bump height (meters)344 r_bump = 50000._wp ! bump radius (meters) 345 h_bump = 2700._wp ! bump height (meters) 332 346 h_oce = gdepw_0(jpk) ! background ocean depth (meters) 333 347 IF(lwp) WRITE(numout,*) ' bump characteristics: ' … … 350 364 idta(:,:) = jpkm1 351 365 DO jk = 1, jpkm1 352 DO jj = 1, jpjdta 353 DO ji = 1, jpidta 354 IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk 355 END DO 356 END DO 366 WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) ) idta(:,:) = jk 357 367 END DO 358 368 ENDIF … … 361 371 ! ! Caution : idta on the global domain: use of jperio, not nperio 362 372 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 363 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1. e0364 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0. e0373 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp 374 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 365 375 ELSEIF( jperio == 2 ) THEN 366 376 idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 ) 367 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0. e0368 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0. e0369 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0. e0377 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 378 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp 379 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp 370 380 ELSE 371 ih = 0 ; zh = 0. e0381 ih = 0 ; zh = 0._wp 372 382 IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce 373 383 idta( : , 1 ) = ih ; zdta( : , 1 ) = zh … … 379 389 ! ! local domain level and meter bathymetries (mbathy,bathy) 380 390 mbathy(:,:) = 0 ! set to zero extra halo points 381 bathy (:,:) = 0. e0! (require for mpp case)391 bathy (:,:) = 0._wp ! (require for mpp case) 382 392 DO jj = 1, nlcj ! interior values 383 393 DO ji = 1, nlci … … 392 402 ! 393 403 IF( ln_zco ) THEN ! zco : read level bathymetry 394 CALL iom_open ( 'bathy_level.nc', inum )395 CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )396 CALL iom_close (inum)404 CALL iom_open ( 'bathy_level.nc', inum ) 405 CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy ) 406 CALL iom_close( inum ) 397 407 mbathy(:,:) = INT( bathy(:,:) ) 398 408 ! ! ===================== 399 409 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 400 410 ! ! ===================== 401 IF( n_cla == 0 ) THEN 402 ! 411 IF( nn_cla == 0 ) THEN 403 412 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 404 413 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) … … 409 418 END DO 410 419 IF(lwp) WRITE(numout,*) 411 IF(lwp) WRITE(numout,*) ' 420 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 412 421 ! 413 422 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open … … 419 428 END DO 420 429 IF(lwp) WRITE(numout,*) 421 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 422 ! 430 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 423 431 ENDIF 424 432 ! … … 427 435 ENDIF 428 436 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 429 CALL iom_open ( 'bathy_meter.nc', inum )430 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )431 CALL iom_close (inum)437 CALL iom_open ( 'bathy_meter.nc', inum ) 438 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 439 CALL iom_close( inum ) 432 440 ! ! ===================== 433 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 434 ! ! ===================== 435 IF( n_cla == 0 ) THEN 436 ! 437 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 438 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 441 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 442 ii0 = 142 ; ii1 = 142 ! ===================== 443 ij0 = 51 ; ij1 = 53 444 DO ji = mi0(ii0), mi1(ii1) ! Close Halmera Strait 445 DO jj = mj0(ij0), mj1(ij1) 446 bathy(ji,jj) = 0._wp 447 END DO 448 END DO 449 IF(lwp) WRITE(numout,*) 450 IF(lwp) WRITE(numout,*) ' orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1 451 ENDIF 452 ! ! ===================== 453 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 454 ! ! ===================== 455 IF( nn_cla == 0 ) THEN 456 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 457 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 439 458 DO ji = mi0(ii0), mi1(ii1) 440 459 DO jj = mj0(ij0), mj1(ij1) 441 bathy(ji,jj) = 284. e0460 bathy(ji,jj) = 284._wp 442 461 END DO 443 462 END DO 444 463 IF(lwp) WRITE(numout,*) 445 IF(lwp) WRITE(numout,*) ' 464 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 446 465 ! 447 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open448 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)466 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 467 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 449 468 DO ji = mi0(ii0), mi1(ii1) 450 469 DO jj = mj0(ij0), mj1(ij1) 451 bathy(ji,jj) = 137. e0470 bathy(ji,jj) = 137._wp 452 471 END DO 453 472 END DO 454 473 IF(lwp) WRITE(numout,*) 455 474 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 456 !457 475 ENDIF 458 476 ! … … 473 491 DO ji = ncsi1(jl), ncsi2(jl) 474 492 mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry 475 bathy (ji,jj) = 0. e0493 bathy (ji,jj) = 0._wp 476 494 END DO 477 495 END DO 478 496 END DO 479 497 ENDIF 480 #if defined key_orca_lev10 481 ! ! ================= ! 482 mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 ! 483 ! ! ================= ! 484 IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' ) 485 #endif 486 ! 498 ! 499 ! ! =========================== ! 500 ! ! set a minimum depth ! 501 ! ! =========================== ! 502 IF( rn_hmin < 0._wp ) THEN ; ik = INT( rn_hmin ) + 1 ! from a nb of level 503 ELSE ; ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 504 ENDIF 505 zhmin = gdepw_0(ik+1) ! minimum depth = ik+1 w-levels 506 WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands 507 ELSE WHERE ; bathy(:,:) = MIN( zhmin , bathy(:,:) ) ! min=zhmin over the oceans 508 END WHERE 509 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 510 ! 487 511 END SUBROUTINE zgr_bat 488 512 … … 601 625 IF( lk_mpp ) THEN 602 626 zbathy(:,:) = FLOAT( mbathy(:,:) ) 603 CALL lbc_lnk( zbathy, 'T', 1. )627 CALL lbc_lnk( zbathy, 'T', 1._wp ) 604 628 mbathy(:,:) = INT( zbathy(:,:) ) 605 629 ENDIF … … 641 665 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab 642 666 zbathy(:,:) = FLOAT( mbathy(:,:) ) 643 CALL lbc_lnk( zbathy, 'T', 1. )667 CALL lbc_lnk( zbathy, 'T', 1._wp ) 644 668 mbathy(:,:) = INT( zbathy(:,:) ) 645 669 ENDIF … … 672 696 673 697 698 SUBROUTINE zgr_bot_level 699 !!---------------------------------------------------------------------- 700 !! *** ROUTINE zgr_bot_level *** 701 !! 702 !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) 703 !! 704 !! ** Method : computes from mbathy with a minimum value of 1 over land 705 !! 706 !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest 707 !! ocean level at t-, u- & v-points 708 !! (min value = 1 over land) 709 !!---------------------------------------------------------------------- 710 INTEGER :: ji, jj ! dummy loop indices 711 REAL(wp), DIMENSION(jpi,jpj) :: zmbk ! 2D workspace 712 !!---------------------------------------------------------------------- 713 ! 714 IF(lwp) WRITE(numout,*) 715 IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 716 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 717 ! 718 mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) 719 ! ! bottom k-index of W-level = mbkt+1 720 DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level 721 DO ji = 1, jpim1 722 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) 723 mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) 724 END DO 725 END DO 726 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 727 zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 728 zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 729 ! 730 END SUBROUTINE zgr_bot_level 731 732 674 733 SUBROUTINE zgr_zco 675 734 !!---------------------------------------------------------------------- … … 678 737 !! ** Purpose : define the z-coordinate system 679 738 !! 680 !! ** Method : set 3D coord. arrays to reference 1D array if lk_zco=T739 !! ** Method : set 3D coord. arrays to reference 1D array 681 740 !!---------------------------------------------------------------------- 682 741 INTEGER :: jk 683 742 !!---------------------------------------------------------------------- 684 743 ! 685 IF( .NOT.lk_zco ) THEN 686 DO jk = 1, jpk 687 fsdept(:,:,jk) = gdept_0(jk) 688 fsdepw(:,:,jk) = gdepw_0(jk) 689 fsde3w(:,:,jk) = gdepw_0(jk) 690 fse3t (:,:,jk) = e3t_0(jk) 691 fse3u (:,:,jk) = e3t_0(jk) 692 fse3v (:,:,jk) = e3t_0(jk) 693 fse3f (:,:,jk) = e3t_0(jk) 694 fse3w (:,:,jk) = e3w_0(jk) 695 fse3uw(:,:,jk) = e3w_0(jk) 696 fse3vw(:,:,jk) = e3w_0(jk) 697 END DO 698 ENDIF 744 DO jk = 1, jpk 745 fsdept(:,:,jk) = gdept_0(jk) 746 fsdepw(:,:,jk) = gdepw_0(jk) 747 fsde3w(:,:,jk) = gdepw_0(jk) 748 fse3t (:,:,jk) = e3t_0(jk) 749 fse3u (:,:,jk) = e3t_0(jk) 750 fse3v (:,:,jk) = e3t_0(jk) 751 fse3f (:,:,jk) = e3t_0(jk) 752 fse3w (:,:,jk) = e3w_0(jk) 753 fse3uw(:,:,jk) = e3w_0(jk) 754 fse3vw(:,:,jk) = e3w_0(jk) 755 END DO 699 756 ! 700 757 END SUBROUTINE zgr_zco 701 758 702 #if defined key_zco703 !!----------------------------------------------------------------------704 !! 'key_zco' : "pure" zco (gdep & e3 are 1D arrays)705 !!----------------------------------------------------------------------706 SUBROUTINE zgr_zps ! Empty routine707 END SUBROUTINE zgr_zps708 SUBROUTINE zgr_sco ! Empty routine709 END SUBROUTINE zgr_sco710 711 #else712 !!----------------------------------------------------------------------713 !! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays)714 !!----------------------------------------------------------------------715 759 716 760 SUBROUTINE zgr_zps … … 764 808 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 765 809 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 766 REAL(wp) :: zmax , zmin ! Maximum and minimum depth810 REAL(wp) :: zmax ! Maximum depth 767 811 REAL(wp) :: zdiff ! temporary scalar 812 REAL(wp) :: zrefdep ! temporary scalar 768 813 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprt ! 3D workspace 769 814 !!--------------------------------------------------------------------- … … 774 819 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 775 820 776 ll_print =.FALSE.! Local variable for debugging777 !! ll_print =.TRUE.821 ll_print = .FALSE. ! Local variable for debugging 822 !! ll_print = .TRUE. 778 823 779 824 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth … … 786 831 ! bathymetry in level (from bathy_meter) 787 832 ! =================== 788 zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) 789 zmin = gdepw_0(4) ! minimum depth = 3 levels 790 791 mbathy(:,:) = jpkm1 ! initialize mbathy to the maximum ocean level available 792 793 ! ! storage of land and island's number (zera and negative values) in mbathy 794 WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) ) 795 796 ! bounded value of bathy 797 !!gm bathy(:,:) = MIN( zmax, MAX( bathy(:,:), zmin ) ) !!gm : simpler as zmin is > 0 798 DO jj = 1, jpj 799 DO ji= 1, jpi 800 IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0 801 ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) ) 802 ENDIF 803 END DO 804 END DO 833 zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) 834 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 835 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 836 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 837 END WHERE 805 838 806 839 ! Compute mbathy for ocean points (i.e. the number of ocean levels) … … 810 843 DO jk = jpkm1, 1, -1 811 844 zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) 812 DO jj = 1, jpj 813 DO ji = 1, jpi 814 IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1 815 END DO 816 END DO 845 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 817 846 END DO 818 847 … … 824 853 e3w (:,:,jk) = e3w_0 (jk) 825 854 END DO 826 hdept(:,:) = gdept(:,:,2 )827 hdepw(:,:) = gdepw(:,:,3 )828 855 ! 829 856 DO jj = 1, jpj … … 835 862 zdepwp = bathy(ji,jj) 836 863 ze3tp = bathy(ji,jj) - gdepw_0(ik) 837 ze3wp = 0.5 * e3w_0(ik) * ( 1.+ ( ze3tp/e3t_0(ik) ) )864 ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) ) 838 865 e3t(ji,jj,ik ) = ze3tp 839 866 e3t(ji,jj,ik+1) = ze3tp … … 855 882 e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 856 883 & / ( gdepw_0( ik+1) - gdepw_0(ik) ) 857 e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &858 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )884 e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) ) & 885 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 859 886 ! ... on ik+1 860 887 e3w (ji,jj,ik+1) = e3t (ji,jj,ik) … … 871 898 ik = mbathy(ji,jj) 872 899 IF( ik > 0 ) THEN ! ocean point only 873 hdept(ji,jj) = gdept(ji,jj,ik )874 hdepw(ji,jj) = gdepw(ji,jj,ik+1)875 900 e3tp (ji,jj) = e3t(ji,jj,ik ) 876 901 e3wp (ji,jj) = e3w(ji,jj,ik ) 877 902 ! test 878 903 zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) 879 IF( zdiff <= 0. .AND. lwp ) THEN904 IF( zdiff <= 0._wp .AND. lwp ) THEN 880 905 it = it + 1 881 906 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj … … 890 915 ! Scale factors and depth at U-, V-, UW and VW-points 891 916 DO jk = 1, jpk ! initialisation to z-scale factors 892 e3u (:,:,jk) 893 e3v (:,:,jk) 894 e3uw(:,:,jk) 895 e3vw(:,:,jk) 917 e3u (:,:,jk) = e3t_0(jk) 918 e3v (:,:,jk) = e3t_0(jk) 919 e3uw(:,:,jk) = e3w_0(jk) 920 e3vw(:,:,jk) = e3w_0(jk) 896 921 END DO 897 922 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 898 923 DO jj = 1, jpjm1 899 924 DO ji = 1, fs_jpim1 ! vector opt. 900 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )901 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )925 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) ) 926 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) ) 902 927 e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 903 928 e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) … … 905 930 END DO 906 931 END DO 907 ! ! Boundary conditions 908 CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. ) 909 CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. ) 932 CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions 933 CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp ) 910 934 ! 911 935 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 912 WHERE( e3u (:,:,jk) == 0. e0) e3u (:,:,jk) = e3t_0(jk)913 WHERE( e3v (:,:,jk) == 0. e0) e3v (:,:,jk) = e3t_0(jk)914 WHERE( e3uw(:,:,jk) == 0. e0) e3uw(:,:,jk) = e3w_0(jk)915 WHERE( e3vw(:,:,jk) == 0. e0) e3vw(:,:,jk) = e3w_0(jk)936 WHERE( e3u (:,:,jk) == 0._wp ) e3u (:,:,jk) = e3t_0(jk) 937 WHERE( e3v (:,:,jk) == 0._wp ) e3v (:,:,jk) = e3t_0(jk) 938 WHERE( e3uw(:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_0(jk) 939 WHERE( e3vw(:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_0(jk) 916 940 END DO 917 941 918 942 ! Scale factor at F-point 919 943 DO jk = 1, jpk ! initialisation to z-scale factors 920 e3f 944 e3f(:,:,jk) = e3t_0(jk) 921 945 END DO 922 946 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors … … 927 951 END DO 928 952 END DO 929 CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions953 CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions 930 954 ! 931 955 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 932 WHERE( e3f(:,:,jk) == 0. e0) e3f(:,:,jk) = e3t_0(jk)956 WHERE( e3f(:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_0(jk) 933 957 END DO 934 958 !!gm bug ? : must be a do loop with mj0,mj1 … … 941 965 942 966 ! Control of the sign 943 IF( MINVAL( e3t (:,:,:) ) <= 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )944 IF( MINVAL( e3w (:,:,:) ) <= 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )945 IF( MINVAL( gdept(:,:,:) ) < 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )946 IF( MINVAL( gdepw(:,:,:) ) < 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )967 IF( MINVAL( e3t (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' ) 968 IF( MINVAL( e3w (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' ) 969 IF( MINVAL( gdept(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) 970 IF( MINVAL( gdepw(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) 947 971 948 972 ! Compute gdep3w (vertical sum of e3w) 949 gdep3w(:,:,1) = 0.5 * e3w(:,:,1)973 gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1) 950 974 DO jk = 2, jpk 951 975 gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) 952 976 END DO 953 977 954 978 ! ! ================= ! 955 979 IF(lwp .AND. ll_print) THEN ! Control print ! … … 979 1003 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 980 1004 ENDIF 981 ! 1005 ! 982 1006 END SUBROUTINE zgr_zps 983 1007 … … 993 1017 !! T-points at integer values (between 1 and jpk) 994 1018 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 995 !! 996 !! Reference : ??? 997 !!---------------------------------------------------------------------- 998 REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate 999 REAL(wp) :: pf ! sigma value 1000 !!---------------------------------------------------------------------- 1001 ! 1002 pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb ) ) & 1019 !!---------------------------------------------------------------------- 1020 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 1021 REAL(wp) :: pf ! sigma value 1022 !!---------------------------------------------------------------------- 1023 ! 1024 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 1003 1025 & - TANH( rn_thetb * rn_theta ) ) & 1004 & * ( COSH( rn_theta ) &1005 & + COSH( rn_theta * ( 2. e0 * rn_thetb - 1.e0 ) ) )&1006 & / ( 2. e0* SINH( rn_theta ) )1026 & * ( COSH( rn_theta ) & 1027 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 1028 & / ( 2._wp * SINH( rn_theta ) ) 1007 1029 ! 1008 1030 END FUNCTION fssig … … 1019 1041 !! T-points at integer values (between 1 and jpk) 1020 1042 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1021 !! 1022 !! Reference : ??? 1023 !!---------------------------------------------------------------------- 1024 REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate 1025 REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient 1026 REAL(wp) :: pf1 ! sigma value 1043 !!---------------------------------------------------------------------- 1044 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 1045 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 1046 REAL(wp) :: pf1 ! sigma value 1027 1047 !!---------------------------------------------------------------------- 1028 1048 ! 1029 1049 IF ( rn_theta == 0 ) then ! uniform sigma 1030 pf1 = - (pk1-0.5) / REAL( jpkm1 )1050 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 1031 1051 ELSE ! stretched sigma 1032 pf1 = ( 1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) +&1033 & pbb * ( (tanh( rn_theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*rn_theta) ) /&1034 & (2*tanh(0.5*rn_theta) ))1052 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 1053 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 1054 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 1035 1055 ENDIF 1036 1056 ! … … 1109 1129 ENDIF 1110 1130 1111 gsigw3 = 0. 0d0 ; gsigt3 = 0.0d0 ; gsi3w3 = 0.0d01112 esigt3 = 0. 0d0 ; esigw3 = 0.0d01113 esigtu3 = 0. 0d0 ; esigtv3 = 0.0d0 ; esigtf3 = 0.0d01114 esigwu3 = 0. 0d0 ; esigwv3 = 0.0d01131 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1132 esigt3 = 0._wp ; esigw3 = 0._wp 1133 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1134 esigwu3 = 0._wp ; esigwv3 = 0._wp 1115 1135 1116 1136 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate … … 1124 1144 DO jj = 1, jpj 1125 1145 DO ji = 1, jpi 1126 IF (bathy(ji,jj) .gt. 0.0) THEN1146 IF( bathy(ji,jj) > 0._wp ) THEN 1127 1147 bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1128 1148 ENDIF … … 1140 1160 ! 1141 1161 ! Smooth the bathymetry (if required) 1142 scosrf(:,:) = 0. e0! ocean surface depth (here zero: no under ice-shelf sea)1162 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1143 1163 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1144 1164 ! 1145 1165 jl = 0 1146 zrmax = 1. e01166 zrmax = 1._wp 1147 1167 ! ! ================ ! 1148 DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax )! Iterative loop !1168 DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1149 1169 ! ! ================ ! 1150 1170 jl = jl + 1 1151 zrmax = 0. e01152 zmsk(:,:) = 0. e01171 zrmax = 0._wp 1172 zmsk(:,:) = 0._wp 1153 1173 DO jj = 1, nlcj 1154 1174 DO ji = 1, nlci … … 1158 1178 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1159 1179 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1160 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1. 01161 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1. 01162 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1. 01163 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1. 01180 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1181 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1182 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1183 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1164 1184 END DO 1165 1185 END DO 1166 1186 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1167 1187 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 1168 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )1188 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp ) 1169 1189 DO jj = 1, nlcj 1170 1190 DO ji = 1, nlci … … 1181 1201 iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) 1182 1202 ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) 1183 IF( zmsk(ji,jj) == 1. 0) THEN1203 IF( zmsk(ji,jj) == 1._wp ) THEN 1184 1204 ztmp(ji,jj) = ( & 1185 1205 & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & 1186 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2. e0+ zenv(iip1,jj )*zmsk(iip1,jj ) &1206 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) & 1187 1207 & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & 1188 1208 & ) / ( & 1189 1209 & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & 1190 & + zmsk(iim1,jj ) + 2. e0+ zmsk(iip1,jj ) &1210 & + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) & 1191 1211 & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & 1192 1212 & ) … … 1197 1217 DO jj = 1, nlcj 1198 1218 DO ji = 1, nlci 1199 IF( zmsk(ji,jj) == 1. 0) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )1219 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1200 1220 END DO 1201 1221 END DO 1202 1222 ! 1203 1223 ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero 1204 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )1224 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1205 1225 DO jj = 1, nlcj 1206 1226 DO ji = 1, nlci 1207 IF( zenv(ji,jj) == 0. e0) zenv(ji,jj) = ztmp(ji,jj)1227 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1208 1228 END DO 1209 1229 END DO … … 1214 1234 ! ! envelop bathymetry saved in hbatt 1215 1235 hbatt(:,:) = zenv(:,:) 1216 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. e0) THEN1236 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1217 1237 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 1218 1238 DO jj = 1, jpj 1219 1239 DO ji = 1, jpi 1220 ztaper = EXP( -(gphit(ji,jj)/8 )**2 )1221 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1.-ztaper)1240 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 ) 1241 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 1222 1242 END DO 1223 1243 END DO … … 1228 1248 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 1229 1249 WRITE(numout,*) 1230 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0. , numout )1250 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout ) 1231 1251 IF( nprint == 1 ) THEN 1232 1252 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) … … 1247 1267 DO jj = 1, jpjm1 1248 1268 DO ji = 1, jpim1 ! NO vector opt. 1249 hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )1250 hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )1251 hbatf(ji,jj) = 0.25 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &1252 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )1269 hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 1270 hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 1271 hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 1272 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 1253 1273 END DO 1254 1274 END DO … … 1256 1276 ! Apply lateral boundary condition 1257 1277 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 1258 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )1278 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp ) 1259 1279 DO jj = 1, jpj 1260 1280 DO ji = 1, jpi 1261 IF( hbatu(ji,jj) == 0. e0) THEN1262 IF( zhbat(ji,jj) == 0. e0) hbatu(ji,jj) = rn_sbot_min1263 IF( zhbat(ji,jj) /= 0. e0) hbatu(ji,jj) = zhbat(ji,jj)1281 IF( hbatu(ji,jj) == 0._wp ) THEN 1282 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 1283 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) 1264 1284 ENDIF 1265 1285 END DO 1266 1286 END DO 1267 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )1287 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp ) 1268 1288 DO jj = 1, jpj 1269 1289 DO ji = 1, jpi 1270 IF( hbatv(ji,jj) == 0. e0) THEN1271 IF( zhbat(ji,jj) == 0. e0) hbatv(ji,jj) = rn_sbot_min1272 IF( zhbat(ji,jj) /= 0. e0) hbatv(ji,jj) = zhbat(ji,jj)1290 IF( hbatv(ji,jj) == 0._wp ) THEN 1291 IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min 1292 IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj) 1273 1293 ENDIF 1274 1294 END DO 1275 1295 END DO 1276 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )1296 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp ) 1277 1297 DO jj = 1, jpj 1278 1298 DO ji = 1, jpi 1279 IF( hbatf(ji,jj) == 0. e0) THEN1280 IF( zhbat(ji,jj) == 0. e0) hbatf(ji,jj) = rn_sbot_min1281 IF( zhbat(ji,jj) /= 0. e0) hbatf(ji,jj) = zhbat(ji,jj)1299 IF( hbatf(ji,jj) == 0._wp ) THEN 1300 IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min 1301 IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj) 1282 1302 ENDIF 1283 1303 END DO … … 1313 1333 DO jj = 1, jpj 1314 1334 1315 IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma 1335 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1336 DO jk = 1, jpk 1337 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1338 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1339 END DO 1340 ELSE ! shallow water, uniform sigma 1341 DO jk = 1, jpk 1342 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1343 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1344 END DO 1345 ENDIF 1346 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1347 ! 1348 DO jk = 1, jpkm1 1349 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1350 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1351 END DO 1352 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) 1353 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 1354 ! 1355 ! Coefficients for vertical depth as the sum of e3w scale factors 1356 gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 1357 DO jk = 2, jpk 1358 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1359 END DO 1360 ! 1316 1361 DO jk = 1, jpk 1317 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1318 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1362 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1363 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1364 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1365 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1366 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1319 1367 END DO 1320 ELSE ! shallow water, uniform sigma 1368 ! 1369 END DO ! for all jj's 1370 END DO ! for all ji's 1371 1372 DO ji = 1, jpi 1373 DO jj = 1, jpj 1321 1374 DO jk = 1, jpk 1322 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1323 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 1375 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & 1376 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1377 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & 1378 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1379 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & 1380 & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & 1381 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1382 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & 1383 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1384 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & 1385 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1386 ! 1387 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1388 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1389 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1390 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1391 ! 1392 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1393 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1394 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1324 1395 END DO 1325 ENDIF 1326 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1327 1328 1329 DO jk = 1, jpkm1 1330 esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1331 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1332 END DO 1333 esigw3(ji,jj,1 ) = 2.0_wp * (gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 )) 1334 esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk)) 1335 1336 ! Coefficients for vertical depth as the sum of e3w scale factors 1337 gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 1338 DO jk = 2, jpk 1339 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1340 END DO 1341 1342 DO jk = 1, jpk 1343 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp) 1344 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp) 1345 gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft) 1346 gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw) 1347 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft) 1348 END DO 1349 1350 ENDDO ! for all jj's 1351 ENDDO ! for all ji's 1352 1353 1354 DO ji=1,jpi 1355 DO jj=1,jpj 1356 1357 DO jk = 1, jpk 1358 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / & 1359 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1360 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / & 1361 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1362 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) + & 1363 hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / & 1364 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1365 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / & 1366 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1367 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / & 1368 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1369 1370 e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1371 e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1372 e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1373 e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1374 ! 1375 e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1376 e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1377 e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1378 END DO 1379 1380 ENDDO 1381 ENDDO 1382 1396 END DO 1397 END DO 1398 ! 1383 1399 ELSE ! not ln_s_sigma 1384 1400 ! 1385 1401 DO jk = 1, jpk 1386 1402 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1387 gsigt(jk) = -fssig( REAL(jk,wp) )1388 END DO 1389 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)1403 gsigt(jk) = -fssig( REAL(jk,wp) ) 1404 END DO 1405 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1390 1406 ! 1391 ! Coefficients for vertical scale factors at w-, t- levels1407 ! Coefficients for vertical scale factors at w-, t- levels 1392 1408 !!gm bug : define it from analytical function, not like juste bellow.... 1393 1409 !!gm or betteroffer the 2 possibilities.... … … 1396 1412 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1397 1413 END DO 1398 esigw( 1 ) = 2. 0_wp * (gsigt(1) - gsigw(1))1399 esigt(jpk) = 2. 0_wp * (gsigt(jpk) - gsigw(jpk))1414 esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) 1415 esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 1400 1416 1401 1417 !!gm original form … … 1407 1423 ! 1408 1424 ! Coefficients for vertical depth as the sum of e3w scale factors 1409 gsi3w(1) = 0.5 * esigw(1)1425 gsi3w(1) = 0.5_wp * esigw(1) 1410 1426 DO jk = 2, jpk 1411 1427 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) … … 1413 1429 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1414 1430 DO jk = 1, jpk 1415 zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)1416 zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)1417 gdept (:,:,jk) = ( scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)1418 gdepw (:,:,jk) = ( scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)1419 gdep3w(:,:,jk) = ( scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoeft)1431 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1432 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1433 gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigt(jk) + hift(:,:)*zcoeft ) 1434 gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsigw(jk) + hift(:,:)*zcoefw ) 1435 gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*gsi3w(jk) + hift(:,:)*zcoeft ) 1420 1436 END DO 1421 1437 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 1423 1439 DO ji = 1, jpi 1424 1440 DO jk = 1, jpk 1425 e3t(ji,jj,jk) =((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))1426 e3u(ji,jj,jk) =((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))1427 e3v(ji,jj,jk) =((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))1428 e3f(ji,jj,jk) =((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))1441 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1442 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1443 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1444 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1429 1445 ! 1430 e3w (ji,jj,jk) =((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))1431 e3uw(ji,jj,jk) =((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))1432 e3vw(ji,jj,jk) =((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))1446 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1447 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1448 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1433 1449 END DO 1434 1450 END DO 1435 1451 END DO 1436 1452 ! 1437 1453 ENDIF ! ln_s_sigma 1438 1454 … … 1457 1473 DO jk = 1, jpkm1 1458 1474 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1459 IF( scobot(ji,jj) == 0. e0) mbathy(ji,jj) = 01475 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 1460 1476 END DO 1461 1477 END DO … … 1463 1479 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 1464 1480 & ' MAX ', MAXVAL( mbathy(:,:) ) 1465 1466 1481 1467 1482 ! ! ============= … … 1523 1538 DO jj = 1, jpj 1524 1539 DO ji = 1, jpi 1525 IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0.) THEN1540 IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1526 1541 WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1527 1542 CALL ctl_stop( ctmp1 ) 1528 1543 ENDIF 1529 IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0.) THEN1544 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1530 1545 WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1531 1546 CALL ctl_stop( ctmp1 ) … … 1538 1553 END SUBROUTINE zgr_sco 1539 1554 1540 #endif1541 1555 1542 1556 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
