- Timestamp:
- 2010-11-25T20:03:49+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r2435 r2436 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 … … 17 17 18 18 !!---------------------------------------------------------------------- 19 !! dom_zgr 20 !! zgr_bat 21 !! zgr_bat_zoom 22 !! zgr_bat_ctl 23 !! zgr_z 24 !! zgr_zco 25 !! zgr_zps 26 !! zgr_sco 27 !! fssig 28 !! dfssig 19 !! dom_zgr : defined the ocean vertical coordinate system 20 !! zgr_bat : bathymetry fields (levels and meters) 21 !! zgr_bat_zoom: modify the bathymetry field if zoom domain 22 !! zgr_bat_ctl : check the bathymetry files 23 !! zgr_z : reference z-coordinate 24 !! zgr_zco : z-coordinate 25 !! zgr_zps : z-coordinate with partial steps 26 !! zgr_sco : s-coordinate 27 !! fssig : sigma coordinate non-dimensional function 28 !! dfssig : derivative of the sigma coordinate function !!gm (currently missing!) 29 29 !!--------------------------------------------------------------------- 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 30 USE oce ! ocean variables 31 USE dom_oce ! ocean domain 32 USE closea ! closed seas 33 USE c1d ! 1D vertical configuration 34 USE in_out_manager ! I/O manager 35 USE iom ! I/O library 36 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 37 USE lib_mpp ! distributed memory computing library 37 38 38 39 IMPLICIT NONE 39 40 PRIVATE 40 41 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 42 PUBLIC dom_zgr ! called by dom_init.F90 43 44 ! !!* Namelist namzgr_sco * 45 REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m) 46 REAL(wp) :: rn_sbot_max = 5250._wp ! maximum depth of s-bottom surface (= ocean depth) (>0) (m) 47 REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20) 48 REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1) 49 REAL(wp) :: rn_rmax = 0.15_wp ! 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.80_wp ! 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._wp ! Critical depth for s-sigma coordinates 54 54 55 55 !! * Substitutions … … 59 59 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 60 60 !! $Id$ 61 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)61 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 62 62 !!---------------------------------------------------------------------- 63 64 63 CONTAINS 65 64 … … 82 81 !!---------------------------------------------------------------------- 83 82 INTEGER :: ioptio = 0 ! temporary integer 84 ! !83 ! 85 84 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco 86 85 !!---------------------------------------------------------------------- 87 86 88 REWIND ( numnam )! Read Namelist namzgr : vertical coordinate'89 READ 87 REWIND( numnam ) ! Read Namelist namzgr : vertical coordinate' 88 READ ( numnam, namzgr ) 90 89 91 90 IF(lwp) THEN ! Control print … … 103 102 IF( ln_zps ) ioptio = ioptio + 1 104 103 IF( ln_sco ) ioptio = ioptio + 1 105 IF 104 IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 106 105 107 106 ! Build the vertical coordinate system … … 113 112 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 114 113 115 !!bug gm control print:116 114 IF( nprint == 1 .AND. lwp ) THEN 117 115 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) … … 130 128 & ' w ', MAXVAL( fse3w(:,:,:) ) 131 129 ENDIF 132 !!!bug gm 133 130 ! 134 131 END SUBROUTINE dom_zgr 135 132 … … 172 169 ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed 173 170 ! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 174 !175 171 IF( ppa1 == pp_to_be_computed .AND. & 176 172 & ppa0 == pp_to_be_computed .AND. & … … 191 187 WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates' 192 188 WRITE(numout,*) ' ~~~~~~~' 193 IF( ppkth == 0. ) THEN189 IF( ppkth == 0._wp ) THEN 194 190 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 195 191 WRITE(numout,*) ' Total depth :', zhmax 196 192 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 197 193 ELSE 198 IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0.) THEN194 IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN 199 195 WRITE(numout,*) ' zsur, za0, za1 computed from ' 200 196 WRITE(numout,*) ' zdzmin = ', zdzmin … … 218 214 ! Reference z-coordinate (depth - scale factor at T- and W-points) 219 215 ! ====================== 220 IF( ppkth == 0. e0) THEN ! uniform vertical grid216 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 221 217 za1 = zhmax / FLOAT(jpk-1) 222 218 DO jk = 1, jpk 223 219 zw = FLOAT( jk ) 224 zt = FLOAT( jk ) + 0.5 220 zt = FLOAT( jk ) + 0.5_wp 225 221 gdepw_0(jk) = ( zw - 1 ) * za1 226 222 gdept_0(jk) = ( zt - 1 ) * za1 … … 232 228 DO jk = 1, jpk 233 229 zw = FLOAT( jk ) 234 zt = FLOAT( jk ) + 0.5 230 zt = FLOAT( jk ) + 0.5_wp 235 231 gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) ) 236 232 gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) ) … … 241 237 DO jk = 1, jpk 242 238 zw = FLOAT( jk ) 243 zt = FLOAT( jk ) + 0.5 239 zt = FLOAT( jk ) + 0.5_wp 244 240 ! Double tanh function 245 241 gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) & … … 253 249 END DO 254 250 ENDIF 255 gdepw_0(1) = 0. e0! force first w-level to be exactly at zero251 gdepw_0(1) = 0._wp ! force first w-level to be exactly at zero 256 252 ENDIF 257 253 258 254 !!gm BUG in s-coordinate this does not work! 259 255 ! deepest/shallowest W level Above/Below ~10m 260 zrefdep = 10. - ( 0.1*MINVAL(e3w_0) )! ref. depth with tolerance (10% of minimum layer thickness)256 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_0 ) ! ref. depth with tolerance (10% of minimum layer thickness) 261 257 nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 262 258 nla10 = nlb10 - 1 ! deepest W level Above ~10m … … 270 266 ENDIF 271 267 DO jk = 1, jpk ! control positivity 272 IF( e3w_0 (jk) <= 0. e0 .OR. e3t_0 (jk) <= 0.e0 ) CALL ctl_stop( 'e3w or e3t =< 0 ' )273 IF( gdepw_0(jk) < 0. e0 .OR. gdept_0(jk) < 0.e0 ) CALL ctl_stop( 'gdepw or gdept < 0 ' )268 IF( e3w_0 (jk) <= 0._wp .OR. e3t_0 (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 ' ) 269 IF( gdepw_0(jk) < 0._wp .OR. gdept_0(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' ) 274 270 END DO 275 271 ! … … 310 306 !! - bathy : meter bathymetry (in meters) 311 307 !!---------------------------------------------------------------------- 312 USE iom313 !!314 308 INTEGER :: ji, jj, jl, jk ! dummy loop indices 315 309 INTEGER :: inum ! temporary logical unit … … 342 336 ii_bump = jpidta / 2 ! i-index of the bump center 343 337 ij_bump = jpjdta / 2 ! j-index of the bump center 344 r_bump = 50000. e0! bump radius (meters)345 h_bump = 2700.e0! bump height (meters)338 r_bump = 50000._wp ! bump radius (meters) 339 h_bump = 2700._wp ! bump height (meters) 346 340 h_oce = gdepw_0(jpk) ! background ocean depth (meters) 347 341 IF(lwp) WRITE(numout,*) ' bump characteristics: ' … … 364 358 idta(:,:) = jpkm1 365 359 DO jk = 1, jpkm1 366 DO jj = 1, jpjdta 367 DO ji = 1, jpidta 368 IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk 369 END DO 370 END DO 360 WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) ) idta(:,:) = jk 371 361 END DO 372 362 ENDIF … … 375 365 ! ! Caution : idta on the global domain: use of jperio, not nperio 376 366 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 377 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1. e0378 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0. e0367 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp 368 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 379 369 ELSEIF( jperio == 2 ) THEN 380 370 idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 ) 381 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0. e0382 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0. e0383 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0. e0371 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 372 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp 373 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp 384 374 ELSE 385 ih = 0 ; zh = 0. e0375 ih = 0 ; zh = 0._wp 386 376 IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce 387 377 idta( : , 1 ) = ih ; zdta( : , 1 ) = zh … … 393 383 ! ! local domain level and meter bathymetries (mbathy,bathy) 394 384 mbathy(:,:) = 0 ! set to zero extra halo points 395 bathy (:,:) = 0. e0! (require for mpp case)385 bathy (:,:) = 0._wp ! (require for mpp case) 396 386 DO jj = 1, nlcj ! interior values 397 387 DO ji = 1, nlci … … 406 396 ! 407 397 IF( ln_zco ) THEN ! zco : read level bathymetry 408 CALL iom_open ( 'bathy_level.nc', inum )409 CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )410 CALL iom_close (inum)398 CALL iom_open ( 'bathy_level.nc', inum ) 399 CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy ) 400 CALL iom_close( inum ) 411 401 mbathy(:,:) = INT( bathy(:,:) ) 412 402 ! ! ===================== … … 439 429 ENDIF 440 430 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 441 CALL iom_open ( 'bathy_meter.nc', inum )442 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )443 CALL iom_close (inum)431 CALL iom_open ( 'bathy_meter.nc', inum ) 432 CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy ) 433 CALL iom_close( inum ) 444 434 ! ! ===================== 445 435 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration … … 448 438 DO ji = mi0(ii0), mi1(ii1) ! Close Halmera Strait 449 439 DO jj = mj0(ij0), mj1(ij1) 450 bathy(ji,jj) = 0. 0440 bathy(ji,jj) = 0._wp 451 441 END DO 452 442 END DO … … 462 452 DO ji = mi0(ii0), mi1(ii1) 463 453 DO jj = mj0(ij0), mj1(ij1) 464 bathy(ji,jj) = 284. e0454 bathy(ji,jj) = 284._wp 465 455 END DO 466 456 END DO … … 472 462 DO ji = mi0(ii0), mi1(ii1) 473 463 DO jj = mj0(ij0), mj1(ij1) 474 bathy(ji,jj) = 137. e0464 bathy(ji,jj) = 137._wp 475 465 END DO 476 466 END DO … … 495 485 DO ji = ncsi1(jl), ncsi2(jl) 496 486 mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry 497 bathy (ji,jj) = 0. e0487 bathy (ji,jj) = 0._wp 498 488 END DO 499 489 END DO 500 490 END DO 501 491 ENDIF 502 #if defined key_orca_lev10503 ! ! ================= !504 mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 !505 ! ! ================= !506 IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )507 #endif508 492 ! ! =============== ! 509 IF( lzoom ) CALL zgr_bat_zoom! Zoom domain !493 IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain ! 510 494 ! ! =============== ! 511 495 ! 512 496 ! ! =================== ! 513 IF( .NOT. lk_c1d ) CALL zgr_bat_ctl! Bathymetry check !497 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check ! 514 498 ! ! =================== ! 515 499 END SUBROUTINE zgr_bat … … 629 613 IF( lk_mpp ) THEN 630 614 zbathy(:,:) = FLOAT( mbathy(:,:) ) 631 CALL lbc_lnk( zbathy, 'T', 1. )615 CALL lbc_lnk( zbathy, 'T', 1._wp ) 632 616 mbathy(:,:) = INT( zbathy(:,:) ) 633 617 ENDIF … … 669 653 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab 670 654 zbathy(:,:) = FLOAT( mbathy(:,:) ) 671 CALL lbc_lnk( zbathy, 'T', 1. )655 CALL lbc_lnk( zbathy, 'T', 1._wp ) 672 656 mbathy(:,:) = INT( zbathy(:,:) ) 673 657 ENDIF … … 792 776 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 793 777 794 ll_print =.FALSE.! Local variable for debugging795 !! ll_print =.TRUE.778 ll_print = .FALSE. ! Local variable for debugging 779 !! ll_print = .TRUE. 796 780 797 781 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth … … 806 790 zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) ) 807 791 808 zrefdep = 25. 792 zrefdep = 25._wp 809 793 nlbelow = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below zrefdep 810 794 IF(lwp) write(numout,*) 'Minimum depth level selected: ', nlbelow … … 814 798 815 799 ! ! storage of land and island's number (zera and negative values) in mbathy 816 WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) )800 WHERE( bathy(:,:) <= 0._wp ) mbathy(:,:) = INT( bathy(:,:) ) 817 801 818 802 ! bounded value of bathy … … 820 804 DO jj = 1, jpj 821 805 DO ji= 1, jpi 822 IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0823 ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) )806 IF( bathy(ji,jj) <= 0._wp ) THEN ; bathy(ji,jj) = 0._wp 807 ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) ) 824 808 ENDIF 825 809 END DO … … 832 816 DO jk = jpkm1, 1, -1 833 817 zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) 834 DO jj = 1, jpj 835 DO ji = 1, jpi 836 IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1 837 END DO 838 END DO 818 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 839 819 END DO 840 820 … … 857 837 zdepwp = bathy(ji,jj) 858 838 ze3tp = bathy(ji,jj) - gdepw_0(ik) 859 ze3wp = 0.5 * e3w_0(ik) * ( 1.+ ( ze3tp/e3t_0(ik) ) )839 ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) ) 860 840 e3t(ji,jj,ik ) = ze3tp 861 841 e3t(ji,jj,ik+1) = ze3tp … … 877 857 e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 878 858 & / ( gdepw_0( ik+1) - gdepw_0(ik) ) 879 e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &880 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )859 e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) ) & 860 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 881 861 ! ... on ik+1 882 862 e3w (ji,jj,ik+1) = e3t (ji,jj,ik) … … 899 879 ! test 900 880 zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) 901 IF( zdiff <= 0. .AND. lwp ) THEN881 IF( zdiff <= 0._wp .AND. lwp ) THEN 902 882 it = it + 1 903 883 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj … … 920 900 DO jj = 1, jpjm1 921 901 DO ji = 1, fs_jpim1 ! vector opt. 922 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )923 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )902 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) ) 903 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) ) 924 904 e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 925 905 e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) … … 927 907 END DO 928 908 END DO 929 ! ! Boundary conditions 930 CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. ) 931 CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. ) 909 CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions 910 CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp ) 932 911 ! 933 912 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 934 WHERE( e3u (:,:,jk) == 0. e0) e3u (:,:,jk) = e3t_0(jk)935 WHERE( e3v (:,:,jk) == 0. e0) e3v (:,:,jk) = e3t_0(jk)936 WHERE( e3uw(:,:,jk) == 0. e0) e3uw(:,:,jk) = e3w_0(jk)937 WHERE( e3vw(:,:,jk) == 0. e0) e3vw(:,:,jk) = e3w_0(jk)913 WHERE( e3u (:,:,jk) == 0._wp ) e3u (:,:,jk) = e3t_0(jk) 914 WHERE( e3v (:,:,jk) == 0._wp ) e3v (:,:,jk) = e3t_0(jk) 915 WHERE( e3uw(:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_0(jk) 916 WHERE( e3vw(:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_0(jk) 938 917 END DO 939 918 940 919 ! Scale factor at F-point 941 920 DO jk = 1, jpk ! initialisation to z-scale factors 942 e3f 921 e3f(:,:,jk) = e3t_0(jk) 943 922 END DO 944 923 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors … … 949 928 END DO 950 929 END DO 951 CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions930 CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions 952 931 ! 953 932 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 954 WHERE( e3f(:,:,jk) == 0. e0) e3f(:,:,jk) = e3t_0(jk)933 WHERE( e3f(:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_0(jk) 955 934 END DO 956 935 !!gm bug ? : must be a do loop with mj0,mj1 … … 963 942 964 943 ! Control of the sign 965 IF( MINVAL( e3t (:,:,:) ) <= 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )966 IF( MINVAL( e3w (:,:,:) ) <= 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )967 IF( MINVAL( gdept(:,:,:) ) < 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )968 IF( MINVAL( gdepw(:,:,:) ) < 0. e0) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )944 IF( MINVAL( e3t (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' ) 945 IF( MINVAL( e3w (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' ) 946 IF( MINVAL( gdept(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) 947 IF( MINVAL( gdepw(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' ) 969 948 970 949 ! Compute gdep3w (vertical sum of e3w) 971 gdep3w(:,:,1) = 0.5 * e3w(:,:,1)950 gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1) 972 951 DO jk = 2, jpk 973 952 gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk) … … 1001 980 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1002 981 ENDIF 1003 982 ! 1004 983 ! ! =============== ! 1005 984 IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain ! 1006 985 ! ! =============== ! 1007 1008 986 ! 1009 987 ! ! =================== ! 1010 988 IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check ! … … 1023 1001 !! T-points at integer values (between 1 and jpk) 1024 1002 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1025 !! 1026 !! Reference : ??? 1027 !!---------------------------------------------------------------------- 1028 REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate 1029 REAL(wp) :: pf ! sigma value 1030 !!---------------------------------------------------------------------- 1031 ! 1032 pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb ) ) & 1003 !!---------------------------------------------------------------------- 1004 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 1005 REAL(wp) :: pf ! sigma value 1006 !!---------------------------------------------------------------------- 1007 ! 1008 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 1033 1009 & - TANH( rn_thetb * rn_theta ) ) & 1034 & * ( COSH( rn_theta ) &1035 & + COSH( rn_theta * ( 2. e0 * rn_thetb - 1.e0 ) ) )&1036 & / ( 2. e0* SINH( rn_theta ) )1010 & * ( COSH( rn_theta ) & 1011 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 1012 & / ( 2._wp * SINH( rn_theta ) ) 1037 1013 ! 1038 1014 END FUNCTION fssig … … 1049 1025 !! T-points at integer values (between 1 and jpk) 1050 1026 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1051 !! 1052 !! Reference : ??? 1053 !!---------------------------------------------------------------------- 1054 REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate 1055 REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient 1056 REAL(wp) :: pf1 ! sigma value 1027 !!---------------------------------------------------------------------- 1028 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 1029 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 1030 REAL(wp) :: pf1 ! sigma value 1057 1031 !!---------------------------------------------------------------------- 1058 1032 ! 1059 1033 IF ( rn_theta == 0 ) then ! uniform sigma 1060 pf1 = - (pk1-0.5) / REAL( jpkm1 )1034 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 1061 1035 ELSE ! stretched sigma 1062 pf1 = ( 1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) +&1063 & pbb * ( (tanh( rn_theta*( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5*rn_theta) ) /&1064 & (2*tanh(0.5*rn_theta) ))1036 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 1037 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 1038 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 1065 1039 ENDIF 1066 1040 ! … … 1139 1113 ENDIF 1140 1114 1141 gsigw3 = 0. 0d0 ; gsigt3 = 0.0d0 ; gsi3w3 = 0.0d01142 esigt3 = 0. 0d0 ; esigw3 = 0.0d01143 esigtu3 = 0. 0d0 ; esigtv3 = 0.0d0 ; esigtf3 = 0.0d01144 esigwu3 = 0. 0d0 ; esigwv3 = 0.0d01115 gsigw3 = 0._wp ; gsigt3 = 0._wp ; gsi3w3 = 0._wp 1116 esigt3 = 0._wp ; esigw3 = 0._wp 1117 esigtu3 = 0._wp ; esigtv3 = 0._wp ; esigtf3 = 0._wp 1118 esigwu3 = 0._wp ; esigwv3 = 0._wp 1145 1119 1146 1120 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate … … 1154 1128 DO jj = 1, jpj 1155 1129 DO ji = 1, jpi 1156 IF (bathy(ji,jj) .gt. 0.0) THEN1130 IF( bathy(ji,jj) > 0._wp ) THEN 1157 1131 bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1158 1132 ENDIF … … 1170 1144 ! 1171 1145 ! Smooth the bathymetry (if required) 1172 scosrf(:,:) = 0. e0! ocean surface depth (here zero: no under ice-shelf sea)1146 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1173 1147 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1174 1148 ! 1175 1149 jl = 0 1176 zrmax = 1. e01150 zrmax = 1._wp 1177 1151 ! ! ================ ! 1178 DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax )! Iterative loop !1152 DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop ! 1179 1153 ! ! ================ ! 1180 1154 jl = jl + 1 1181 zrmax = 0. e01182 zmsk(:,:) = 0. e01155 zrmax = 0._wp 1156 zmsk(:,:) = 0._wp 1183 1157 DO jj = 1, nlcj 1184 1158 DO ji = 1, nlci … … 1188 1162 zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 1189 1163 zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) ) 1190 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1. 01191 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1. 01192 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1. 01193 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1. 01164 IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1165 IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp 1166 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp 1167 IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp 1194 1168 END DO 1195 1169 END DO 1196 1170 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 1197 1171 ! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX) 1198 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )1172 ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp ) 1199 1173 DO jj = 1, nlcj 1200 1174 DO ji = 1, nlci … … 1211 1185 iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci) 1212 1186 ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj) 1213 IF( zmsk(ji,jj) == 1. 0) THEN1187 IF( zmsk(ji,jj) == 1._wp ) THEN 1214 1188 ztmp(ji,jj) = ( & 1215 1189 & zenv(iim1,ijp1)*zmsk(iim1,ijp1) + zenv(ji,ijp1)*zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) & 1216 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2. e0+ zenv(iip1,jj )*zmsk(iip1,jj ) &1190 & + zenv(iim1,jj )*zmsk(iim1,jj ) + zenv(ji,jj )* 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) & 1217 1191 & + zenv(iim1,ijm1)*zmsk(iim1,ijm1) + zenv(ji,ijm1)*zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) & 1218 1192 & ) / ( & 1219 1193 & zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) & 1220 & + zmsk(iim1,jj ) + 2. e0+ zmsk(iip1,jj ) &1194 & + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) & 1221 1195 & + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) & 1222 1196 & ) … … 1227 1201 DO jj = 1, nlcj 1228 1202 DO ji = 1, nlci 1229 IF( zmsk(ji,jj) == 1. 0) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )1203 IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) ) 1230 1204 END DO 1231 1205 END DO 1232 1206 ! 1233 1207 ! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero 1234 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )1208 ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp ) 1235 1209 DO jj = 1, nlcj 1236 1210 DO ji = 1, nlci 1237 IF( zenv(ji,jj) == 0. e0) zenv(ji,jj) = ztmp(ji,jj)1211 IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj) 1238 1212 END DO 1239 1213 END DO … … 1244 1218 ! ! envelop bathymetry saved in hbatt 1245 1219 hbatt(:,:) = zenv(:,:) 1246 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0. e0) THEN1220 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 1247 1221 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 1248 1222 DO jj = 1, jpj 1249 1223 DO ji = 1, jpi 1250 ztaper = EXP( -(gphit(ji,jj)/8 )**2 )1251 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1.-ztaper)1224 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 ) 1225 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 1252 1226 END DO 1253 1227 END DO … … 1258 1232 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters' 1259 1233 WRITE(numout,*) 1260 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0. , numout )1234 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout ) 1261 1235 IF( nprint == 1 ) THEN 1262 1236 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) ) … … 1277 1251 DO jj = 1, jpjm1 1278 1252 DO ji = 1, jpim1 ! NO vector opt. 1279 hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )1280 hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )1281 hbatf(ji,jj) = 0.25 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &1282 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )1253 hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 1254 hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 1255 hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 1256 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 1283 1257 END DO 1284 1258 END DO … … 1286 1260 ! Apply lateral boundary condition 1287 1261 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 1288 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )1262 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp ) 1289 1263 DO jj = 1, jpj 1290 1264 DO ji = 1, jpi 1291 IF( hbatu(ji,jj) == 0. e0) THEN1292 IF( zhbat(ji,jj) == 0. e0) hbatu(ji,jj) = rn_sbot_min1293 IF( zhbat(ji,jj) /= 0. e0) hbatu(ji,jj) = zhbat(ji,jj)1265 IF( hbatu(ji,jj) == 0._wp ) THEN 1266 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 1267 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) 1294 1268 ENDIF 1295 1269 END DO 1296 1270 END DO 1297 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )1271 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp ) 1298 1272 DO jj = 1, jpj 1299 1273 DO ji = 1, jpi 1300 IF( hbatv(ji,jj) == 0. e0) THEN1301 IF( zhbat(ji,jj) == 0. e0) hbatv(ji,jj) = rn_sbot_min1302 IF( zhbat(ji,jj) /= 0. e0) hbatv(ji,jj) = zhbat(ji,jj)1274 IF( hbatv(ji,jj) == 0._wp ) THEN 1275 IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min 1276 IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj) 1303 1277 ENDIF 1304 1278 END DO 1305 1279 END DO 1306 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )1280 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp ) 1307 1281 DO jj = 1, jpj 1308 1282 DO ji = 1, jpi 1309 IF( hbatf(ji,jj) == 0. e0) THEN1310 IF( zhbat(ji,jj) == 0. e0) hbatf(ji,jj) = rn_sbot_min1311 IF( zhbat(ji,jj) /= 0. e0) hbatf(ji,jj) = zhbat(ji,jj)1283 IF( hbatf(ji,jj) == 0._wp ) THEN 1284 IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min 1285 IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj) 1312 1286 ENDIF 1313 1287 END DO … … 1343 1317 DO jj = 1, jpj 1344 1318 1345 IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma 1319 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 1320 DO jk = 1, jpk 1321 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1322 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1323 END DO 1324 ELSE ! shallow water, uniform sigma 1325 DO jk = 1, jpk 1326 gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 1327 gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 1328 END DO 1329 ENDIF 1330 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1331 ! 1332 DO jk = 1, jpkm1 1333 esigt3(ji,jj,jk ) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1334 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1335 END DO 1336 esigw3(ji,jj,1 ) = 2._wp * ( gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ) ) 1337 esigt3(ji,jj,jpk) = 2._wp * ( gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk) ) 1338 ! 1339 ! Coefficients for vertical depth as the sum of e3w scale factors 1340 gsi3w3(ji,jj,1) = 0.5_wp * esigw3(ji,jj,1) 1341 DO jk = 2, jpk 1342 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1343 END DO 1344 ! 1346 1345 DO jk = 1, jpk 1347 gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 1348 gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 1346 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1347 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 1348 gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 1349 gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 1350 gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 1349 1351 END DO 1350 ELSE ! shallow water, uniform sigma 1352 ! 1353 END DO ! for all jj's 1354 END DO ! for all ji's 1355 1356 DO ji = 1, jpi 1357 DO jj = 1, jpj 1351 1358 DO jk = 1, jpk 1352 gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 1353 gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 1359 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) & 1360 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1361 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) & 1362 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1363 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) & 1364 & + hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) & 1365 & / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1366 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) & 1367 & / ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1368 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) & 1369 & / ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1370 ! 1371 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigt3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1372 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1373 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1374 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1375 ! 1376 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*esigw3 (ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1377 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1378 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1) ) 1354 1379 END DO 1355 ENDIF 1356 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk) 1357 1358 1359 DO jk = 1, jpkm1 1360 esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk) 1361 esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk) 1362 END DO 1363 esigw3(ji,jj,1 ) = 2.0_wp * (gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 )) 1364 esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk)) 1365 1366 ! Coefficients for vertical depth as the sum of e3w scale factors 1367 gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1) 1368 DO jk = 2, jpk 1369 gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk) 1370 END DO 1371 1372 DO jk = 1, jpk 1373 zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp) 1374 zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp) 1375 gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigt3(ji,jj,jk)+rn_hc*zcoeft) 1376 gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsigw3(ji,jj,jk)+rn_hc*zcoefw) 1377 gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)*gsi3w3(ji,jj,jk)+rn_hc*zcoeft) 1378 END DO 1379 1380 ENDDO ! for all jj's 1381 ENDDO ! for all ji's 1382 1383 1384 DO ji=1,jpi 1385 DO jj=1,jpj 1386 1387 DO jk = 1, jpk 1388 esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) ) / & 1389 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1390 esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji,jj+1)*esigt3(ji,jj+1,jk) ) / & 1391 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1392 esigtf3(ji,jj,jk) = ( hbatt(ji,jj)*esigt3(ji,jj,jk)+hbatt(ji+1,jj)*esigt3(ji+1,jj,jk) + & 1393 hbatt(ji,jj+1)*esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)*esigt3(ji+1,jj+1,jk) ) / & 1394 ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 1395 esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji+1,jj)*esigw3(ji+1,jj,jk) ) / & 1396 ( hbatt(ji,jj)+hbatt(ji+1,jj) ) 1397 esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*esigw3(ji,jj,jk)+hbatt(ji,jj+1)*esigw3(ji,jj+1,jk) ) / & 1398 ( hbatt(ji,jj)+hbatt(ji,jj+1) ) 1399 1400 e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1401 e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1402 e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1403 e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1404 ! 1405 e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1406 e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1407 e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1)) 1408 END DO 1409 1410 ENDDO 1411 ENDDO 1412 1380 END DO 1381 END DO 1382 ! 1413 1383 ELSE ! not ln_s_sigma 1414 1384 ! 1415 1385 DO jk = 1, jpk 1416 1386 gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 1417 gsigt(jk) = -fssig( REAL(jk,wp) )1418 END DO 1419 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)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) 1420 1390 ! 1421 ! Coefficients for vertical scale factors at w-, t- levels1391 ! Coefficients for vertical scale factors at w-, t- levels 1422 1392 !!gm bug : define it from analytical function, not like juste bellow.... 1423 1393 !!gm or betteroffer the 2 possibilities.... … … 1426 1396 esigw(jk+1) = gsigt(jk+1) - gsigt(jk) 1427 1397 END DO 1428 esigw( 1 ) = 2. 0_wp * (gsigt(1) - gsigw(1))1429 esigt(jpk) = 2. 0_wp * (gsigt(jpk) - gsigw(jpk))1398 esigw( 1 ) = 2._wp * ( gsigt(1 ) - gsigw(1 ) ) 1399 esigt(jpk) = 2._wp * ( gsigt(jpk) - gsigw(jpk) ) 1430 1400 1431 1401 !!gm original form … … 1437 1407 ! 1438 1408 ! Coefficients for vertical depth as the sum of e3w scale factors 1439 gsi3w(1) = 0.5 * esigw(1)1409 gsi3w(1) = 0.5_wp * esigw(1) 1440 1410 DO jk = 2, jpk 1441 1411 gsi3w(jk) = gsi3w(jk-1) + esigw(jk) … … 1443 1413 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 1444 1414 DO jk = 1, jpk 1445 zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)1446 zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)1447 gdept (:,:,jk) = ( scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigt(jk)+hift(:,:)*zcoeft)1448 gdepw (:,:,jk) = ( scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsigw(jk)+hift(:,:)*zcoefw)1449 gdep3w(:,:,jk) = ( scosrf(:,:)+(hbatt(:,:)-hift(:,:))*gsi3w(jk)+hift(:,:)*zcoeft)1415 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 1416 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 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 ) 1450 1420 END DO 1451 1421 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 1453 1423 DO ji = 1, jpi 1454 1424 DO jk = 1, jpk 1455 e3t(ji,jj,jk) =((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))1456 e3u(ji,jj,jk) =((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))1457 e3v(ji,jj,jk) =((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))1458 e3f(ji,jj,jk) =((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))1425 e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1426 e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1427 e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1428 e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 1459 1429 ! 1460 e3w (ji,jj,jk) =((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))1461 e3uw(ji,jj,jk) =((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))1462 e3vw(ji,jj,jk) =((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))1430 e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 1431 e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 1432 e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 1463 1433 END DO 1464 1434 END DO 1465 1435 END DO 1466 1436 ! 1467 1437 ENDIF ! ln_s_sigma 1468 1438 … … 1487 1457 DO jk = 1, jpkm1 1488 1458 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 1489 IF( scobot(ji,jj) == 0. e0) mbathy(ji,jj) = 01459 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 1490 1460 END DO 1491 1461 END DO … … 1498 1468 IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain 1499 1469 ! ! =========== 1500 1501 1470 ! 1502 1471 ! ! =================== ! 1503 1472 IF( .NOT. lk_c1d ) CALL zgr_bat_ctl ! Bathymetry check ! 1504 1473 ! ! =================== ! 1505 1474 ! 1506 1475 ! ! ============= 1507 1476 IF(lwp) THEN ! Control print … … 1562 1531 DO jj = 1, jpj 1563 1532 DO ji = 1, jpi 1564 IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0.) THEN1533 IF( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN 1565 1534 WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1566 1535 CALL ctl_stop( ctmp1 ) 1567 1536 ENDIF 1568 IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0.) THEN1537 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN 1569 1538 WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1570 1539 CALL ctl_stop( ctmp1 )
Note: See TracChangeset
for help on using the changeset viewer.
