Changeset 7646 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r6492 r7646 22 22 23 23 !!---------------------------------------------------------------------- 24 !! dom_zgr : defined the ocean vertical coordinate system 25 !! zgr_bat : bathymetry fields (levels and meters) 26 !! zgr_bat_zoom : modify the bathymetry field if zoom domain 27 !! zgr_bat_ctl : check the bathymetry files 28 !! zgr_bot_level: deepest ocean level for t-, u, and v-points 29 !! zgr_z : reference z-coordinate 30 !! zgr_zco : z-coordinate 31 !! zgr_zps : z-coordinate with partial steps 32 !! zgr_sco : s-coordinate 33 !! fssig : tanh stretch function 34 !! fssig1 : Song and Haidvogel 1994 stretch function 35 !! fgamma : Siddorn and Furner 2012 stretching function 24 !! dom_zgr : read or set the ocean vertical coordinate system 25 !! zgr_read : read the vertical information in the domain configuration file 26 !! zgr_top_bot : ocean top and bottom level for t-, u, and v-points with 1 as minimum value 36 27 !!--------------------------------------------------------------------- 37 USE oce 38 USE dom_oce 39 USE wet_dry ! wetting and drying40 USE closea ! closed seas41 USE c1d ! 1D vertical configuration28 USE oce ! ocean variables 29 USE dom_oce ! ocean domain 30 USE usrdef_zgr ! user defined vertical coordinate system 31 USE depth_e3 ! depth <=> e3 32 USE wet_dry, ONLY: ln_wd, ht_wd 42 33 ! 43 USE in_out_manager 44 USE iom 45 USE lbclnk 46 USE lib_mpp 47 USE wrk_nemo 48 USE timing 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 38 USE wrk_nemo ! Memory allocation 39 USE timing ! Timing 49 40 50 41 IMPLICIT NONE … … 52 43 53 44 PUBLIC dom_zgr ! called by dom_init.F90 54 55 ! !!* Namelist namzgr_sco *56 LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)57 LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)58 !59 REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)60 REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)61 REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)62 REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates63 ! Song and Haidvogel 1994 stretching parameters64 REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)65 REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)66 REAL(wp) :: rn_bb ! stretching parameter67 ! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)68 ! Siddorn and Furner stretching parameters69 LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch70 REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)71 REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord72 REAL(wp) :: rn_zs ! depth of surface grid box73 ! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b74 REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb75 REAL(wp) :: rn_zb_b ! offset for calculating Zb76 45 77 46 !! * Substitutions … … 84 53 CONTAINS 85 54 86 SUBROUTINE dom_zgr 55 SUBROUTINE dom_zgr( k_top, k_bot ) 87 56 !!---------------------------------------------------------------------- 88 57 !! *** ROUTINE dom_zgr *** … … 101 70 !! ** Action : define gdep., e3., mbathy and bathy 102 71 !!---------------------------------------------------------------------- 103 INTEGER :: ioptio, ibat ! local integer 104 INTEGER :: ios 105 ! 106 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 72 INTEGER, DIMENSION(:,:), INTENT(out) :: k_top, k_bot ! ocean first and last level indices 73 ! 74 INTEGER :: jk ! dummy loop index 75 INTEGER :: ioptio, ibat, ios ! local integer 76 REAL(wp) :: zrefdep ! depth of the reference level (~10m) 107 77 !!---------------------------------------------------------------------- 108 78 ! 109 79 IF( nn_timing == 1 ) CALL timing_start('dom_zgr') 110 80 ! 111 REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate112 READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )113 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )114 115 REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate116 READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )117 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )118 IF(lwm) WRITE ( numond, namzgr )119 120 81 IF(lwp) THEN ! Control print 121 82 WRITE(numout,*) 122 83 WRITE(numout,*) 'dom_zgr : vertical coordinate' 123 84 WRITE(numout,*) '~~~~~~~' 124 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 85 ENDIF 86 87 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 88 89 90 IF( ln_read_cfg ) THEN !== read in mesh_mask.nc file ==! 91 IF(lwp) WRITE(numout,*) 92 IF(lwp) WRITE(numout,*) ' Read vertical mesh in ', TRIM( cn_domcfg ), ' file' 93 ! 94 CALL zgr_read ( ln_zco , ln_zps , ln_sco, ln_isfcav, & 95 & gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth 96 & gdept_0 , gdepw_0 , & ! gridpoints depth 97 & e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors 98 & e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors 99 & k_top , k_bot ) ! 1st & last ocean level 100 ! 101 ELSE !== User defined configuration ==! 102 IF(lwp) WRITE(numout,*) 103 IF(lwp) WRITE(numout,*) ' User defined vertical mesh (usr_def_zgr)' 104 ! 105 CALL usr_def_zgr( ln_zco , ln_zps , ln_sco, ln_isfcav, & 106 & gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth 107 & gdept_0 , gdepw_0 , & ! gridpoints depth 108 & e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors 109 & e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors 110 & k_top , k_bot ) ! 1st & last ocean level 111 ! 112 ENDIF 113 ! 114 !!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears 115 ! Compute gde3w_0 (vertical sum of e3w) 116 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 117 DO jk = 2, jpk 118 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 119 END DO 120 ! 121 IF(lwp) THEN ! Control print 122 WRITE(numout,*) 123 WRITE(numout,*) ' Type of vertical coordinate (read in ', TRIM( cn_domcfg ), ' file or set in userdef_nam) :' 125 124 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 126 125 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 127 126 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 128 127 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 129 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 130 ENDIF 131 132 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 128 ENDIF 133 129 134 130 ioptio = 0 ! Check Vertical coordinate options … … 137 133 IF( ln_sco ) ioptio = ioptio + 1 138 134 IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' ) 139 ! 140 ioptio = 0 141 IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1 142 IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1 143 IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' ) 144 ! 145 ! Build the vertical coordinate system 146 ! ------------------------------------ 147 CALL zgr_z ! Reference z-coordinate system (always called) 148 CALL zgr_bat ! Bathymetry fields (levels and meters) 149 IF( lk_c1d ) CALL lbc_lnk( bathy , 'T', 1._wp ) ! 1D config.: same bathy value over the 3x3 domain 150 IF( ln_zco ) CALL zgr_zco ! z-coordinate 151 IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate 152 IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate 153 ! 154 ! final adjustment of mbathy & check 155 ! ----------------------------------- 156 IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain 157 IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points 158 CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points 159 CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points 160 ! 161 IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain 162 ibat = mbathy(2,2) 163 mbathy(:,:) = ibat 164 END IF 135 136 137 ! ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top) 138 CALL zgr_top_bot( k_top, k_bot ) ! with a minimum value set to 1 139 140 141 ! ! deepest/shallowest W level Above/Below ~10m 142 !!gm BUG in s-coordinate this does not work! 143 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness) 144 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 145 nla10 = nlb10 - 1 ! deepest W level Above ~10m 146 !!gm end bug 165 147 ! 166 148 IF( nprint == 1 .AND. lwp ) THEN 167 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 149 WRITE(numout,*) ' MIN val k_top ', MINVAL( k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) ) 150 WRITE(numout,*) ' MIN val k_bot ', MINVAL( k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) ) 168 151 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 169 152 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) … … 186 169 187 170 188 SUBROUTINE zgr_z 189 !!---------------------------------------------------------------------- 190 !! *** ROUTINE zgr_z *** 191 !! 192 !! ** Purpose : set the depth of model levels and the resulting 193 !! vertical scale factors. 194 !! 195 !! ** Method : z-coordinate system (use in all type of coordinate) 196 !! The depth of model levels is defined from an analytical 197 !! function the derivative of which gives the scale factors. 198 !! both depth and scale factors only depend on k (1d arrays). 199 !! w-level: gdepw_1d = gdep(k) 200 !! e3w_1d(k) = dk(gdep)(k) = e3(k) 201 !! t-level: gdept_1d = gdep(k+0.5) 202 !! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5) 203 !! 204 !! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m) 205 !! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m) 206 !! 207 !! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766. 208 !!---------------------------------------------------------------------- 209 INTEGER :: jk ! dummy loop indices 210 REAL(wp) :: zt, zw ! temporary scalars 211 REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in 212 REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90 213 REAL(wp) :: zrefdep ! depth of the reference level (~10m) 214 REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters 215 !!---------------------------------------------------------------------- 216 ! 217 IF( nn_timing == 1 ) CALL timing_start('zgr_z') 218 ! 219 ! Set variables from parameters 220 ! ------------------------------ 221 zkth = ppkth ; zacr = ppacr 222 zdzmin = ppdzmin ; zhmax = pphmax 223 zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters 224 225 ! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed 226 ! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr 227 IF( ppa1 == pp_to_be_computed .AND. & 228 & ppa0 == pp_to_be_computed .AND. & 229 & ppsur == pp_to_be_computed ) THEN 230 ! 231 #if defined key_agrif 232 za1 = ( ppdzmin - pphmax / FLOAT(jpkdta-1) ) & 233 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpkdta-1) * ( LOG( COSH( (jpkdta - ppkth) / ppacr) )& 234 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 235 #else 236 za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) & 237 & / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) & 238 & - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) ) 239 #endif 240 za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr ) 241 zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) ) 242 ELSE 243 za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur 244 za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter 245 ENDIF 246 247 IF(lwp) THEN ! Parameter print 171 SUBROUTINE zgr_read( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate 172 & pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate 173 & pdept , pdepw , & ! 3D t & w-points depth 174 & pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors 175 & pe3w , pe3uw , pe3vw , & ! - - - 176 & k_top , k_bot ) ! top & bottom ocean level 177 !!--------------------------------------------------------------------- 178 !! *** ROUTINE zgr_read *** 179 !! 180 !! ** Purpose : Read the vertical information in the domain configuration file 181 !! 182 !!---------------------------------------------------------------------- 183 LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags 184 LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag 185 REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m] 186 REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m] 187 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m] 188 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m] 189 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - - 190 INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level 191 ! 192 INTEGER :: jk ! dummy loop index 193 INTEGER :: inum ! local logical unit 194 REAL(WP) :: z_zco, z_zps, z_sco, z_cav 195 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 196 !!---------------------------------------------------------------------- 197 ! 198 IF(lwp) THEN 248 199 WRITE(numout,*) 249 WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates' 250 WRITE(numout,*) ' ~~~~~~~' 251 IF( ppkth == 0._wp ) THEN 252 WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers' 253 WRITE(numout,*) ' Total depth :', zhmax 254 #if defined key_agrif 255 WRITE(numout,*) ' Layer thickness:', zhmax/(jpkdta-1) 256 #else 257 WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1) 258 #endif 259 ELSE 260 IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN 261 WRITE(numout,*) ' zsur, za0, za1 computed from ' 262 WRITE(numout,*) ' zdzmin = ', zdzmin 263 WRITE(numout,*) ' zhmax = ', zhmax 264 ENDIF 265 WRITE(numout,*) ' Value of coefficients for vertical mesh:' 266 WRITE(numout,*) ' zsur = ', zsur 267 WRITE(numout,*) ' za0 = ', za0 268 WRITE(numout,*) ' za1 = ', za1 269 WRITE(numout,*) ' zkth = ', zkth 270 WRITE(numout,*) ' zacr = ', zacr 271 IF( ldbletanh ) THEN 272 WRITE(numout,*) ' (Double tanh za2 = ', za2 273 WRITE(numout,*) ' parameters) zkth2= ', zkth2 274 WRITE(numout,*) ' zacr2= ', zacr2 275 ENDIF 200 WRITE(numout,*) ' zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file' 201 WRITE(numout,*) ' ~~~~~~~~' 202 ENDIF 203 ! 204 CALL iom_open( cn_domcfg, inum ) 205 ! 206 ! !* type of vertical coordinate 207 CALL iom_get( inum, 'ln_zco' , z_zco ) 208 CALL iom_get( inum, 'ln_zps' , z_zps ) 209 CALL iom_get( inum, 'ln_sco' , z_sco ) 210 IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF 211 IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF 212 IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF 213 ! 214 ! !* ocean cavities under iceshelves 215 CALL iom_get( inum, 'ln_isfcav', z_cav ) 216 IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF 217 ! 218 ! !* vertical scale factors 219 CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate 220 CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d ) 221 ! 222 CALL iom_get( inum, jpdom_data, 'e3t_0' , pe3t , lrowattr=ln_use_jattr ) ! 3D coordinate 223 CALL iom_get( inum, jpdom_data, 'e3u_0' , pe3u , lrowattr=ln_use_jattr ) 224 CALL iom_get( inum, jpdom_data, 'e3v_0' , pe3v , lrowattr=ln_use_jattr ) 225 CALL iom_get( inum, jpdom_data, 'e3f_0' , pe3f , lrowattr=ln_use_jattr ) 226 CALL iom_get( inum, jpdom_data, 'e3w_0' , pe3w , lrowattr=ln_use_jattr ) 227 CALL iom_get( inum, jpdom_data, 'e3uw_0' , pe3uw , lrowattr=ln_use_jattr ) 228 CALL iom_get( inum, jpdom_data, 'e3vw_0' , pe3vw , lrowattr=ln_use_jattr ) 229 ! 230 ! !* depths 231 ! !- old depth definition (obsolescent feature) 232 IF( iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0 .AND. & 233 & iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0 .AND. & 234 & iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0 .AND. & 235 & iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0 ) THEN 236 CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', & 237 & ' depths at t- and w-points read in the domain configuration file') 238 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d ) 239 CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d ) 240 CALL iom_get( inum, jpdom_data , 'gdept_0' , pdept , lrowattr=ln_use_jattr ) 241 CALL iom_get( inum, jpdom_data , 'gdepw_0' , pdepw , lrowattr=ln_use_jattr ) 242 ! 243 ELSE !- depths computed from e3. scale factors 244 CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) ! 1D reference depth 245 CALL e3_to_depth( pe3t , pe3w , pdept , pdepw ) ! 3D depths 246 IF(lwp) THEN 247 WRITE(numout,*) 248 WRITE(numout,*) ' Reference 1D z-coordinate depth and scale factors:' 249 WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" ) 250 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk ) 276 251 ENDIF 277 252 ENDIF 278 279 280 ! Reference z-coordinate (depth - scale factor at T- and W-points) 281 ! ====================== 282 IF( ppkth == 0._wp ) THEN ! uniform vertical grid 283 #if defined key_agrif 284 za1 = zhmax / FLOAT(jpkdta-1) 285 #else 286 za1 = zhmax / FLOAT(jpk-1) 287 #endif 288 DO jk = 1, jpk 289 zw = FLOAT( jk ) 290 zt = FLOAT( jk ) + 0.5_wp 291 gdepw_1d(jk) = ( zw - 1 ) * za1 292 gdept_1d(jk) = ( zt - 1 ) * za1 293 e3w_1d (jk) = za1 294 e3t_1d (jk) = za1 295 END DO 296 ELSE ! Madec & Imbard 1996 function 297 IF( .NOT. ldbletanh ) THEN 298 DO jk = 1, jpk 299 zw = REAL( jk , wp ) 300 zt = REAL( jk , wp ) + 0.5_wp 301 gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) ) 302 gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) ) 303 e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr ) 304 e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr ) 305 END DO 306 ELSE 307 DO jk = 1, jpk 308 zw = FLOAT( jk ) 309 zt = FLOAT( jk ) + 0.5_wp 310 ! Double tanh function 311 gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) & 312 & + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) ) 313 gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) & 314 & + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) ) 315 e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) & 316 & + za2 * TANH( (zw-zkth2) / zacr2 ) 317 e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) & 318 & + za2 * TANH( (zt-zkth2) / zacr2 ) 319 END DO 320 ENDIF 321 gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero 322 ENDIF 323 324 IF ( ln_isfcav ) THEN 325 ! need to be like this to compute the pressure gradient with ISF. If not, level beneath the ISF are not aligned (sum(e3t) /= depth) 326 ! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively 327 DO jk = 1, jpkm1 328 e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk) 329 END DO 330 e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO 331 332 DO jk = 2, jpk 333 e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1) 334 END DO 335 e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1)) 336 END IF 337 338 !!gm BUG in s-coordinate this does not work! 339 ! deepest/shallowest W level Above/Below ~10m 340 zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness) 341 nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m 342 nla10 = nlb10 - 1 ! deepest W level Above ~10m 343 !!gm end bug 344 345 IF(lwp) THEN ! control print 346 WRITE(numout,*) 347 WRITE(numout,*) ' Reference z-coordinate depth and scale factors:' 348 WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" ) 349 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk ) 350 ENDIF 351 DO jk = 1, jpk ! control positivity 352 IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' ) 353 IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' ) 354 END DO 355 ! 356 IF( nn_timing == 1 ) CALL timing_stop('zgr_z') 357 ! 358 END SUBROUTINE zgr_z 359 360 361 SUBROUTINE zgr_bat 362 !!---------------------------------------------------------------------- 363 !! *** ROUTINE zgr_bat *** 364 !! 365 !! ** Purpose : set bathymetry both in levels and meters 366 !! 367 !! ** Method : read or define mbathy and bathy arrays 368 !! * level bathymetry: 369 !! The ocean basin geometry is given by a two-dimensional array, 370 !! mbathy, which is defined as follow : 371 !! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level 372 !! at t-point (ji,jj). 373 !! = 0 over the continental t-point. 374 !! The array mbathy is checked to verified its consistency with 375 !! model option. in particular: 376 !! mbathy must have at least 1 land grid-points (mbathy<=0) 377 !! along closed boundary. 378 !! mbathy must be cyclic IF jperio=1. 379 !! mbathy must be lower or equal to jpk-1. 380 !! isolated ocean grid points are suppressed from mbathy 381 !! since they are only connected to remaining 382 !! ocean through vertical diffusion. 383 !! ntopo=-1 : rectangular channel or bassin with a bump 384 !! ntopo= 0 : flat rectangular channel or basin 385 !! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file 386 !! bathy is read in 'bathy_meter.nc' NetCDF file 387 !! 388 !! ** Action : - mbathy: level bathymetry (in level index) 389 !! - bathy : meter bathymetry (in meters) 390 !!---------------------------------------------------------------------- 391 INTEGER :: ji, jj, jk ! dummy loop indices 392 INTEGER :: inum ! temporary logical unit 393 INTEGER :: ierror ! error flag 394 INTEGER :: ii_bump, ij_bump, ih ! bump center position 395 INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices 396 REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics 397 REAL(wp) :: zi, zj, zh, zhmin ! local scalars 398 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data 399 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data 400 !!---------------------------------------------------------------------- 401 ! 402 IF( nn_timing == 1 ) CALL timing_start('zgr_bat') 403 ! 404 IF(lwp) WRITE(numout,*) 405 IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry' 406 IF(lwp) WRITE(numout,*) ' ~~~~~~~' 407 ! ! ================== ! 408 IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand ! 409 ! ! ================== ! 410 ! ! global domain level and meter bathymetry (idta,zdta) 411 ! 412 ALLOCATE( idta(jpidta,jpjdta), STAT=ierror ) 413 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' ) 414 ALLOCATE( zdta(jpidta,jpjdta), STAT=ierror ) 415 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' ) 416 ! 417 IF( ntopo == 0 ) THEN ! flat basin 418 IF(lwp) WRITE(numout,*) 419 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin' 420 IF( rn_bathy > 0.01 ) THEN 421 IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist' 422 zdta(:,:) = rn_bathy 423 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 424 idta(:,:) = jpkm1 425 ELSE ! z-coordinate (zco or zps): step-like topography 426 idta(:,:) = jpkm1 427 DO jk = 1, jpkm1 428 WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk 429 END DO 430 ENDIF 431 ELSE 432 IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)' 433 idta(:,:) = jpkm1 ! before last level 434 zdta(:,:) = gdepw_1d(jpk) ! last w-point depth 435 h_oce = gdepw_1d(jpk) 436 ENDIF 437 ELSE ! bump centered in the basin 438 IF(lwp) WRITE(numout,*) 439 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump' 440 ii_bump = jpidta / 2 ! i-index of the bump center 441 ij_bump = jpjdta / 2 ! j-index of the bump center 442 r_bump = 50000._wp ! bump radius (meters) 443 h_bump = 2700._wp ! bump height (meters) 444 h_oce = gdepw_1d(jpk) ! background ocean depth (meters) 445 IF(lwp) WRITE(numout,*) ' bump characteristics: ' 446 IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump 447 IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters' 448 IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index' 449 IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters' 450 ! 451 DO jj = 1, jpjdta ! zdta : 452 DO ji = 1, jpidta 453 zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump 454 zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump 455 zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) ) 456 END DO 457 END DO 458 ! ! idta : 459 IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk 460 idta(:,:) = jpkm1 461 ELSE ! z-coordinate (zco or zps): step-like topography 462 idta(:,:) = jpkm1 463 DO jk = 1, jpkm1 464 WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk 465 END DO 466 ENDIF 467 ENDIF 468 ! ! set GLOBAL boundary conditions 469 ! ! Caution : idta on the global domain: use of jperio, not nperio 470 IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN 471 idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp 472 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 473 ELSEIF( jperio == 2 ) THEN 474 idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 ) 475 idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp 476 idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp 477 idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp 478 ELSE 479 ih = 0 ; zh = 0._wp 480 IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce 481 idta( : , 1 ) = ih ; zdta( : , 1 ) = zh 482 idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh 483 idta( 1 , : ) = ih ; zdta( 1 , : ) = zh 484 idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh 485 ENDIF 486 487 ! ! local domain level and meter bathymetries (mbathy,bathy) 488 mbathy(:,:) = 0 ! set to zero extra halo points 489 bathy (:,:) = 0._wp ! (require for mpp case) 490 DO jj = 1, nlcj ! interior values 491 DO ji = 1, nlci 492 mbathy(ji,jj) = idta( mig(ji), mjg(jj) ) 493 bathy (ji,jj) = zdta( mig(ji), mjg(jj) ) 494 END DO 495 END DO 496 risfdep(:,:)=0.e0 497 misfdep(:,:)=1 498 ! 499 DEALLOCATE( idta, zdta ) 500 ! 501 ! ! ================ ! 502 ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain) 503 ! ! ================ ! 504 ! 505 IF( ln_zco ) THEN ! zco : read level bathymetry 506 CALL iom_open ( 'bathy_level.nc', inum ) 507 CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy ) 508 CALL iom_close( inum ) 509 mbathy(:,:) = INT( bathy(:,:) ) 510 ! initialisation isf variables 511 risfdep(:,:)=0._wp ; misfdep(:,:)=1 512 ! ! ===================== 513 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 514 ! ! ===================== 515 ! 516 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 517 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 518 DO ji = mi0(ii0), mi1(ii1) 519 DO jj = mj0(ij0), mj1(ij1) 520 mbathy(ji,jj) = 15 521 END DO 522 END DO 523 IF(lwp) WRITE(numout,*) 524 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 525 ! 526 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 527 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 528 DO ji = mi0(ii0), mi1(ii1) 529 DO jj = mj0(ij0), mj1(ij1) 530 mbathy(ji,jj) = 12 531 END DO 532 END DO 533 IF(lwp) WRITE(numout,*) 534 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 535 ! 536 ENDIF 537 ! 538 ENDIF 539 IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry 540 CALL iom_open ( 'bathy_meter.nc', inum ) 541 IF ( ln_isfcav ) THEN 542 CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. ) 543 ELSE 544 CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) 545 END IF 546 CALL iom_close( inum ) 547 ! 548 ! initialisation isf variables 549 risfdep(:,:)=0._wp ; misfdep(:,:)=1 550 ! 551 IF ( ln_isfcav ) THEN 552 CALL iom_open ( 'isf_draft_meter.nc', inum ) 553 CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep ) 554 CALL iom_close( inum ) 555 WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp 556 557 ! set grounded point to 0 558 ! (a treshold could be set here if needed, or set it offline based on the grounded fraction) 559 WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin ) 560 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 561 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 562 END WHERE 563 END IF 564 ! 565 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 566 ! 567 ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open 568 ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995) 569 DO ji = mi0(ii0), mi1(ii1) 570 DO jj = mj0(ij0), mj1(ij1) 571 bathy(ji,jj) = 284._wp 572 END DO 573 END DO 574 IF(lwp) WRITE(numout,*) 575 IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0 576 ! 577 ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open 578 ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995) 579 DO ji = mi0(ii0), mi1(ii1) 580 DO jj = mj0(ij0), mj1(ij1) 581 bathy(ji,jj) = 137._wp 582 END DO 583 END DO 584 IF(lwp) WRITE(numout,*) 585 IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0 586 ! 587 ENDIF 588 ! 589 ENDIF 590 ! ! =============== ! 591 ELSE ! error ! 592 ! ! =============== ! 593 WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo 594 CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) ) 595 ENDIF 596 ! 597 IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==! 598 ! 599 IF ( .not. ln_sco ) THEN !== set a minimum depth ==! 600 IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level 601 ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth 602 ENDIF 603 zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels 604 WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands 605 ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans 606 END WHERE 607 IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik 608 ENDIF 609 ! 610 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat') 611 ! 612 END SUBROUTINE zgr_bat 613 614 615 SUBROUTINE zgr_bat_zoom 616 !!---------------------------------------------------------------------- 617 !! *** ROUTINE zgr_bat_zoom *** 618 !! 619 !! ** Purpose : - Close zoom domain boundary if necessary 620 !! - Suppress Med Sea from ORCA R2 and R05 arctic zoom 621 !! 622 !! ** Method : 623 !! 624 !! ** Action : - update mbathy: level bathymetry (in level index) 625 !!---------------------------------------------------------------------- 626 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 627 !!---------------------------------------------------------------------- 628 ! 629 IF(lwp) WRITE(numout,*) 630 IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain' 631 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~' 632 ! 633 ! Zoom domain 634 ! =========== 635 ! 636 ! Forced closed boundary if required 637 IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0 638 IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0 639 IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0 640 IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0 641 ! 642 ! Configuration specific domain modifications 643 ! (here, ORCA arctic configuration: suppress Med Sea) 644 IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN 645 SELECT CASE ( jp_cfg ) 646 ! ! ======================= 647 CASE ( 2 ) ! ORCA_R2 configuration 648 ! ! ======================= 649 IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea' 650 ii0 = 141 ; ii1 = 162 ! Sea box i,j indices 651 ij0 = 98 ; ij1 = 110 652 ! ! ======================= 653 CASE ( 05 ) ! ORCA_R05 configuration 654 ! ! ======================= 655 IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea' 656 ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe 657 ij0 = 314 ; ij1 = 370 658 END SELECT 659 ! 660 mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe 661 ! 662 ENDIF 663 ! 664 END SUBROUTINE zgr_bat_zoom 665 666 667 SUBROUTINE zgr_bat_ctl 668 !!---------------------------------------------------------------------- 669 !! *** ROUTINE zgr_bat_ctl *** 670 !! 671 !! ** Purpose : check the bathymetry in levels 672 !! 673 !! ** Method : The array mbathy is checked to verified its consistency 674 !! with the model options. in particular: 675 !! mbathy must have at least 1 land grid-points (mbathy<=0) 676 !! along closed boundary. 677 !! mbathy must be cyclic IF jperio=1. 678 !! mbathy must be lower or equal to jpk-1. 679 !! isolated ocean grid points are suppressed from mbathy 680 !! since they are only connected to remaining 681 !! ocean through vertical diffusion. 682 !! C A U T I O N : mbathy will be modified during the initializa- 683 !! tion phase to become the number of non-zero w-levels of a water 684 !! column, with a minimum value of 1. 685 !! 686 !! ** Action : - update mbathy: level bathymetry (in level index) 687 !! - update bathy : meter bathymetry (in meters) 688 !!---------------------------------------------------------------------- 689 INTEGER :: ji, jj, jl ! dummy loop indices 690 INTEGER :: icompt, ibtest, ikmax ! temporary integers 691 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 692 !!---------------------------------------------------------------------- 693 ! 694 IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl') 695 ! 696 CALL wrk_alloc( jpi, jpj, zbathy ) 697 ! 698 IF(lwp) WRITE(numout,*) 699 IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry' 700 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 701 ! ! Suppress isolated ocean grid points 702 IF(lwp) WRITE(numout,*) 703 IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points' 704 IF(lwp) WRITE(numout,*)' -----------------------------------' 705 icompt = 0 706 DO jl = 1, 2 707 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 708 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 709 mbathy(jpi,:) = mbathy( 2 ,:) 710 ENDIF 711 DO jj = 2, jpjm1 712 DO ji = 2, jpim1 713 ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), & 714 & mbathy(ji,jj-1), mbathy(ji,jj+1) ) 715 IF( ibtest < mbathy(ji,jj) ) THEN 716 IF(lwp) WRITE(numout,*) ' the number of ocean level at ', & 717 & 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest 718 mbathy(ji,jj) = ibtest 719 icompt = icompt + 1 720 ENDIF 721 END DO 722 END DO 723 END DO 724 IF( lk_mpp ) CALL mpp_sum( icompt ) 725 IF( icompt == 0 ) THEN 726 IF(lwp) WRITE(numout,*)' no isolated ocean grid points' 727 ELSE 728 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed' 729 ENDIF 730 IF( lk_mpp ) THEN 731 zbathy(:,:) = FLOAT( mbathy(:,:) ) 732 CALL lbc_lnk( zbathy, 'T', 1._wp ) 733 mbathy(:,:) = INT( zbathy(:,:) ) 734 ENDIF 735 ! ! East-west cyclic boundary conditions 736 IF( nperio == 0 ) THEN 737 IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio 738 IF( lk_mpp ) THEN 739 IF( nbondi == -1 .OR. nbondi == 2 ) THEN 740 IF( jperio /= 1 ) mbathy(1,:) = 0 741 ENDIF 742 IF( nbondi == 1 .OR. nbondi == 2 ) THEN 743 IF( jperio /= 1 ) mbathy(nlci,:) = 0 744 ENDIF 745 ELSE 746 IF( ln_zco .OR. ln_zps ) THEN 747 mbathy( 1 ,:) = 0 748 mbathy(jpi,:) = 0 749 ELSE 750 mbathy( 1 ,:) = jpkm1 751 mbathy(jpi,:) = jpkm1 752 ENDIF 753 ENDIF 754 ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 755 IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio 756 mbathy( 1 ,:) = mbathy(jpim1,:) 757 mbathy(jpi,:) = mbathy( 2 ,:) 758 ELSEIF( nperio == 2 ) THEN 759 IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio 760 ELSE 761 IF(lwp) WRITE(numout,*) ' e r r o r' 762 IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio 763 ! STOP 'dom_mba' 764 ENDIF 765 ! Boundary condition on mbathy 766 IF( .NOT.lk_mpp ) THEN 767 !!gm !!bug ??? think about it ! 768 ! ... mono- or macro-tasking: T-point, >0, 2D array, no slab 769 zbathy(:,:) = FLOAT( mbathy(:,:) ) 770 CALL lbc_lnk( zbathy, 'T', 1._wp ) 771 mbathy(:,:) = INT( zbathy(:,:) ) 772 ENDIF 773 ! Number of ocean level inferior or equal to jpkm1 774 ikmax = 0 775 DO jj = 1, jpj 776 DO ji = 1, jpi 777 ikmax = MAX( ikmax, mbathy(ji,jj) ) 778 END DO 779 END DO 780 !!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ??? 781 IF( ikmax > jpkm1 ) THEN 782 IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1' 783 IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry' 784 ELSE IF( ikmax < jpkm1 ) THEN 785 IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1' 786 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 787 ENDIF 788 ! 789 CALL wrk_dealloc( jpi, jpj, zbathy ) 790 ! 791 IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl') 792 ! 793 END SUBROUTINE zgr_bat_ctl 794 795 796 SUBROUTINE zgr_bot_level 797 !!---------------------------------------------------------------------- 798 !! *** ROUTINE zgr_bot_level *** 253 ! 254 ! !* ocean top and bottom level 255 CALL iom_get( inum, jpdom_data, 'top_level' , z2d , lrowattr=ln_use_jattr ) ! 1st wet T-points (ISF) 256 k_top(:,:) = INT( z2d(:,:) ) 257 CALL iom_get( inum, jpdom_data, 'bottom_level' , z2d , lrowattr=ln_use_jattr ) ! last wet T-points 258 k_bot(:,:) = INT( z2d(:,:) ) 259 ! 260 ! bathymetry with orography (wetting and drying only) 261 IF( ln_wd ) CALL iom_get( inum, jpdom_data, 'ht_wd' , ht_wd , lrowattr=ln_use_jattr ) 262 ! 263 CALL iom_close( inum ) 264 ! 265 END SUBROUTINE zgr_read 266 267 268 SUBROUTINE zgr_top_bot( k_top, k_bot ) 269 !!---------------------------------------------------------------------- 270 !! *** ROUTINE zgr_top_bot *** 799 271 !! 800 272 !! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays) 801 273 !! 802 !! ** Method : computes from mbathy with a minimum value of 1 over land 803 !! 274 !! ** Method : computes from k_top and k_bot with a minimum value of 1 over land 275 !! 276 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 277 !! ocean level at t-, u- & v-points 278 !! (min value = 1) 804 279 !! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest 805 280 !! ocean level at t-, u- & v-points 806 281 !! (min value = 1 over land) 807 282 !!---------------------------------------------------------------------- 283 INTEGER , DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! top & bottom ocean level indices 284 ! 808 285 INTEGER :: ji, jj ! dummy loop indices 809 REAL(wp), POINTER, DIMENSION(:,:) :: z mbk810 !!---------------------------------------------------------------------- 811 ! 812 IF( nn_timing == 1 ) CALL timing_start('zgr_ bot_level')813 ! 814 CALL wrk_alloc( jpi, jpj, zmbk )286 REAL(wp), POINTER, DIMENSION(:,:) :: zk 287 !!---------------------------------------------------------------------- 288 ! 289 IF( nn_timing == 1 ) CALL timing_start('zgr_top_bot') 290 ! 291 CALL wrk_alloc( jpi,jpj, zk ) 815 292 ! 816 293 IF(lwp) WRITE(numout,*) 817 IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels ' 818 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 819 ! 820 mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land) 294 IF(lwp) WRITE(numout,*) ' zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels ' 295 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~' 296 ! 297 mikt(:,:) = MAX( k_top(:,:) , 1 ) ! top ocean k-index of T-level (=1 over land) 298 ! 299 mbkt(:,:) = MAX( k_bot(:,:) , 1 ) ! bottom ocean k-index of T-level (=1 over land) 821 300 822 ! ! bottom k-index of W-level = mbkt+1 823 DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level 301 ! ! N.B. top k-index of W-level = mikt 302 ! ! bottom k-index of W-level = mbkt+1 303 DO jj = 1, jpjm1 824 304 DO ji = 1, jpim1 305 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 306 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 307 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 308 ! 825 309 mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) ) 826 310 mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) ) … … 828 312 END DO 829 313 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 830 zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 831 zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 ) 832 ! 833 CALL wrk_dealloc( jpi, jpj, zmbk ) 834 ! 835 IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level') 836 ! 837 END SUBROUTINE zgr_bot_level 838 839 840 SUBROUTINE zgr_top_level 841 !!---------------------------------------------------------------------- 842 !! *** ROUTINE zgr_top_level *** 843 !! 844 !! ** Purpose : defines the vertical index of ocean top (mik. arrays) 845 !! 846 !! ** Method : computes from misfdep with a minimum value of 1 847 !! 848 !! ** Action : mikt, miku, mikv : vertical indices of the shallowest 849 !! ocean level at t-, u- & v-points 850 !! (min value = 1) 851 !!---------------------------------------------------------------------- 852 INTEGER :: ji, jj ! dummy loop indices 853 REAL(wp), POINTER, DIMENSION(:,:) :: zmik 854 !!---------------------------------------------------------------------- 855 ! 856 IF( nn_timing == 1 ) CALL timing_start('zgr_top_level') 857 ! 858 CALL wrk_alloc( jpi, jpj, zmik ) 859 ! 860 IF(lwp) WRITE(numout,*) 861 IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels ' 862 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~' 863 ! 864 mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1) 865 ! ! top k-index of W-level (=mikt) 866 DO jj = 1, jpjm1 ! top k-index of U- (U-) level 867 DO ji = 1, jpim1 868 miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) ) 869 mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) ) 870 mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) ) 871 END DO 872 END DO 873 874 ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk 875 zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk(zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 ) 876 zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk(zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 ) 877 zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk(zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 ) 878 ! 879 CALL wrk_dealloc( jpi, jpj, zmik ) 880 ! 881 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_level') 882 ! 883 END SUBROUTINE zgr_top_level 884 885 886 SUBROUTINE zgr_zco 887 !!---------------------------------------------------------------------- 888 !! *** ROUTINE zgr_zco *** 889 !! 890 !! ** Purpose : define the reference z-coordinate system 891 !! 892 !! ** Method : set 3D coord. arrays to reference 1D array 893 !!---------------------------------------------------------------------- 894 INTEGER :: jk 895 !!---------------------------------------------------------------------- 896 ! 897 IF( nn_timing == 1 ) CALL timing_start('zgr_zco') 898 ! 899 DO jk = 1, jpk 900 gdept_0(:,:,jk) = gdept_1d(jk) 901 gdepw_0(:,:,jk) = gdepw_1d(jk) 902 gde3w_0(:,:,jk) = gdepw_1d(jk) 903 e3t_0 (:,:,jk) = e3t_1d (jk) 904 e3u_0 (:,:,jk) = e3t_1d (jk) 905 e3v_0 (:,:,jk) = e3t_1d (jk) 906 e3f_0 (:,:,jk) = e3t_1d (jk) 907 e3w_0 (:,:,jk) = e3w_1d (jk) 908 e3uw_0 (:,:,jk) = e3w_1d (jk) 909 e3vw_0 (:,:,jk) = e3w_1d (jk) 910 END DO 911 ! 912 IF( nn_timing == 1 ) CALL timing_stop('zgr_zco') 913 ! 914 END SUBROUTINE zgr_zco 915 916 917 SUBROUTINE zgr_zps 918 !!---------------------------------------------------------------------- 919 !! *** ROUTINE zgr_zps *** 920 !! 921 !! ** Purpose : the depth and vertical scale factor in partial step 922 !! reference z-coordinate case 923 !! 924 !! ** Method : Partial steps : computes the 3D vertical scale factors 925 !! of T-, U-, V-, W-, UW-, VW and F-points that are associated with 926 !! a partial step representation of bottom topography. 927 !! 928 !! The reference depth of model levels is defined from an analytical 929 !! function the derivative of which gives the reference vertical 930 !! scale factors. 931 !! From depth and scale factors reference, we compute there new value 932 !! with partial steps on 3d arrays ( i, j, k ). 933 !! 934 !! w-level: gdepw_0(i,j,k) = gdep(k) 935 !! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k) 936 !! t-level: gdept_0(i,j,k) = gdep(k+0.5) 937 !! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5) 938 !! 939 !! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), 940 !! we find the mbathy index of the depth at each grid point. 941 !! This leads us to three cases: 942 !! 943 !! - bathy = 0 => mbathy = 0 944 !! - 1 < mbathy < jpkm1 945 !! - bathy > gdepw_0(jpk) => mbathy = jpkm1 946 !! 947 !! Then, for each case, we find the new depth at t- and w- levels 948 !! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw- 949 !! and f-points. 950 !! 951 !! This routine is given as an example, it must be modified 952 !! following the user s desiderata. nevertheless, the output as 953 !! well as the way to compute the model levels and scale factors 954 !! must be respected in order to insure second order accuracy 955 !! schemes. 956 !! 957 !! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives 958 !! - - - - - - - gdept_0, gdepw_0 and e3. are positives 959 !! 960 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 961 !!---------------------------------------------------------------------- 962 INTEGER :: ji, jj, jk ! dummy loop indices 963 INTEGER :: ik, it, ikb, ikt ! temporary integers 964 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 965 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 966 REAL(wp) :: zdiff ! temporary scalar 967 REAL(wp) :: zmax ! temporary scalar 968 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 969 !!--------------------------------------------------------------------- 970 ! 971 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 972 ! 973 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 974 ! 975 IF(lwp) WRITE(numout,*) 976 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 977 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 978 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 979 980 ! bathymetry in level (from bathy_meter) 981 ! =================== 982 zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) ) 983 bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat) 984 WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0 985 ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level 986 END WHERE 987 988 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 989 ! find the number of ocean levels such that the last level thickness 990 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where 991 ! e3t_1d is the reference level thickness 992 DO jk = jpkm1, 1, -1 993 zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 994 WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1 995 END DO 996 997 ! Scale factors and depth at T- and W-points 998 DO jk = 1, jpk ! intitialization to the reference z-coordinate 999 gdept_0(:,:,jk) = gdept_1d(jk) 1000 gdepw_0(:,:,jk) = gdepw_1d(jk) 1001 e3t_0 (:,:,jk) = e3t_1d (jk) 1002 e3w_0 (:,:,jk) = e3w_1d (jk) 1003 END DO 1004 1005 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 1006 IF ( ln_isfcav ) CALL zgr_isf 1007 1008 ! Scale factors and depth at T- and W-points 1009 IF ( .NOT. ln_isfcav ) THEN 1010 DO jj = 1, jpj 1011 DO ji = 1, jpi 1012 ik = mbathy(ji,jj) 1013 IF( ik > 0 ) THEN ! ocean point only 1014 ! max ocean level case 1015 IF( ik == jpkm1 ) THEN 1016 zdepwp = bathy(ji,jj) 1017 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1018 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1019 e3t_0(ji,jj,ik ) = ze3tp 1020 e3t_0(ji,jj,ik+1) = ze3tp 1021 e3w_0(ji,jj,ik ) = ze3wp 1022 e3w_0(ji,jj,ik+1) = ze3tp 1023 gdepw_0(ji,jj,ik+1) = zdepwp 1024 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1025 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1026 ! 1027 ELSE ! standard case 1028 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1029 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1030 ENDIF 1031 !gm Bug? check the gdepw_1d 1032 ! ... on ik 1033 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1034 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1035 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1036 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1037 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1038 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1039 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1040 ! ... on ik+1 1041 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1042 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1043 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1044 ENDIF 1045 ENDIF 1046 END DO 1047 END DO 1048 ! 1049 it = 0 1050 DO jj = 1, jpj 1051 DO ji = 1, jpi 1052 ik = mbathy(ji,jj) 1053 IF( ik > 0 ) THEN ! ocean point only 1054 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1055 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1056 ! test 1057 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1058 IF( zdiff <= 0._wp .AND. lwp ) THEN 1059 it = it + 1 1060 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1061 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1062 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1063 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1064 ENDIF 1065 ENDIF 1066 END DO 1067 END DO 1068 END IF 1069 ! 1070 ! Scale factors and depth at U-, V-, UW and VW-points 1071 DO jk = 1, jpk ! initialisation to z-scale factors 1072 e3u_0 (:,:,jk) = e3t_1d(jk) 1073 e3v_0 (:,:,jk) = e3t_1d(jk) 1074 e3uw_0(:,:,jk) = e3w_1d(jk) 1075 e3vw_0(:,:,jk) = e3w_1d(jk) 1076 END DO 1077 1078 DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors 1079 DO jj = 1, jpjm1 1080 DO ji = 1, fs_jpim1 ! vector opt. 1081 e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) ) 1082 e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) ) 1083 e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) ) 1084 e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) ) 1085 END DO 1086 END DO 1087 END DO 1088 IF ( ln_isfcav ) THEN 1089 ! (ISF) define e3uw (adapted for 2 cells in the water column) 1090 DO jj = 2, jpjm1 1091 DO ji = 2, fs_jpim1 ! vector opt. 1092 ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj)) 1093 ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj)) 1094 IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) & 1095 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) ) 1096 ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1)) 1097 ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1)) 1098 IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) & 1099 & - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) ) 1100 END DO 1101 END DO 1102 END IF 1103 1104 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions 1105 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 1106 ! 1107 1108 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1109 WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk) 1110 WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk) 1111 WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk) 1112 WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk) 1113 END DO 1114 1115 ! Scale factor at F-point 1116 DO jk = 1, jpk ! initialisation to z-scale factors 1117 e3f_0(:,:,jk) = e3t_1d(jk) 1118 END DO 1119 DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors 1120 DO jj = 1, jpjm1 1121 DO ji = 1, fs_jpim1 ! vector opt. 1122 e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) ) 1123 END DO 1124 END DO 1125 END DO 1126 CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions 1127 ! 1128 DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries) 1129 WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk) 1130 END DO 1131 !!gm bug ? : must be a do loop with mj0,mj1 1132 ! 1133 e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2 1134 e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:) 1135 e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:) 1136 e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:) 1137 e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:) 1138 1139 ! Control of the sign 1140 IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' ) 1141 IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' ) 1142 IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' ) 1143 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1144 1145 ! Compute gde3w_0 (vertical sum of e3w) 1146 IF ( ln_isfcav ) THEN ! if cavity 1147 WHERE( misfdep == 0 ) misfdep = 1 1148 DO jj = 1,jpj 1149 DO ji = 1,jpi 1150 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1151 DO jk = 2, misfdep(ji,jj) 1152 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1153 END DO 1154 IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1155 DO jk = misfdep(ji,jj) + 1, jpk 1156 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1157 END DO 1158 END DO 1159 END DO 1160 ELSE ! no cavity 1161 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1162 DO jk = 2, jpk 1163 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1164 END DO 1165 END IF 1166 ! 1167 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1168 ! 1169 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1170 ! 1171 END SUBROUTINE zgr_zps 1172 1173 1174 SUBROUTINE zgr_isf 1175 !!---------------------------------------------------------------------- 1176 !! *** ROUTINE zgr_isf *** 1177 !! 1178 !! ** Purpose : check the bathymetry in levels 1179 !! 1180 !! ** Method : THe water column have to contained at least 2 cells 1181 !! Bathymetry and isfdraft are modified (dig/close) to respect 1182 !! this criterion. 1183 !! 1184 !! ** Action : - test compatibility between isfdraft and bathy 1185 !! - bathy and isfdraft are modified 1186 !!---------------------------------------------------------------------- 1187 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1188 INTEGER :: ik, it ! temporary integers 1189 INTEGER :: icompt, ibtest ! (ISF) 1190 INTEGER :: ibtestim1, ibtestip1 ! (ISF) 1191 INTEGER :: ibtestjm1, ibtestjp1 ! (ISF) 1192 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1193 REAL(wp) :: zmax ! Maximum and minimum depth 1194 REAL(wp) :: zbathydiff ! isf temporary scalar 1195 REAL(wp) :: zrisfdepdiff ! isf temporary scalar 1196 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1197 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1198 REAL(wp) :: zdiff ! temporary scalar 1199 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1200 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) 1201 !!--------------------------------------------------------------------- 1202 ! 1203 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1204 ! 1205 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1206 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1207 1208 ! (ISF) compute misfdep 1209 WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1 1210 ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level 1211 END WHERE 1212 1213 ! Compute misfdep for ocean points (i.e. first wet level) 1214 ! find the first ocean level such that the first level thickness 1215 ! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where 1216 ! e3t_0 is the reference level thickness 1217 DO jk = 2, jpkm1 1218 zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1219 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1220 END DO 1221 WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) ) 1222 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1223 END WHERE 1224 1225 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1226 WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1) 1227 misfdep = 0; risfdep = 0.0_wp; 1228 mbathy = 0; bathy = 0.0_wp; 1229 END WHERE 1230 WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp) 1231 misfdep = 0; risfdep = 0.0_wp; 1232 mbathy = 0; bathy = 0.0_wp; 1233 END WHERE 1234 1235 ! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation 1236 icompt = 0 1237 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1238 DO jl = 1, 10 1239 ! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments) 1240 WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin) 1241 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1242 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp 1243 END WHERE 1244 WHERE (mbathy(:,:) <= 0) 1245 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1246 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1247 END WHERE 1248 IF( lk_mpp ) THEN 1249 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1250 CALL lbc_lnk( zbathy, 'T', 1. ) 1251 misfdep(:,:) = INT( zbathy(:,:) ) 1252 1253 CALL lbc_lnk( risfdep,'T', 1. ) 1254 CALL lbc_lnk( bathy, 'T', 1. ) 1255 1256 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1257 CALL lbc_lnk( zbathy, 'T', 1. ) 1258 mbathy(:,:) = INT( zbathy(:,:) ) 1259 ENDIF 1260 IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN 1261 misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west 1262 misfdep(jpi,:) = misfdep( 2 ,:) 1263 mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west 1264 mbathy(jpi,:) = mbathy( 2 ,:) 1265 ENDIF 1266 1267 ! split last cell if possible (only where water column is 2 cell or less) 1268 ! if coupled to ice sheet, we do not modify the bathymetry (can be discuss). 1269 IF ( .NOT. ln_iscpl) THEN 1270 DO jk = jpkm1, 1, -1 1271 zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1272 WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy) 1273 mbathy(:,:) = jk 1274 bathy(:,:) = zmax 1275 END WHERE 1276 END DO 1277 END IF 1278 1279 ! split top cell if possible (only where water column is 2 cell or less) 1280 DO jk = 2, jpkm1 1281 zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat ) 1282 WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy) 1283 misfdep(:,:) = jk 1284 risfdep(:,:) = zmax 1285 END WHERE 1286 END DO 1287 1288 1289 ! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition 1290 DO jj = 1, jpj 1291 DO ji = 1, jpi 1292 ! find the minimum change option: 1293 ! test bathy 1294 IF (risfdep(ji,jj) > 1) THEN 1295 IF ( .NOT. ln_iscpl ) THEN 1296 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1297 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1298 zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) & 1299 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1300 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1301 IF (zbathydiff <= zrisfdepdiff) THEN 1302 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) 1303 mbathy(ji,jj)= mbathy(ji,jj) + 1 1304 ELSE 1305 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1306 misfdep(ji,jj) = misfdep(ji,jj) - 1 1307 END IF 1308 ENDIF 1309 ELSE 1310 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1311 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ) 1312 misfdep(ji,jj) = misfdep(ji,jj) - 1 1313 END IF 1314 END IF 1315 END IF 1316 END DO 1317 END DO 1318 1319 ! At least 2 levels for water thickness at T, U, and V point. 1320 DO jj = 1, jpj 1321 DO ji = 1, jpi 1322 ! find the minimum change option: 1323 ! test bathy 1324 IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1325 IF ( .NOT. ln_iscpl ) THEN 1326 zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1327 & + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) 1328 zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) & 1329 & - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1330 IF (zbathydiff <= zrisfdepdiff) THEN 1331 mbathy(ji,jj) = mbathy(ji,jj) + 1 1332 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1333 ELSE 1334 misfdep(ji,jj)= misfdep(ji,jj) - 1 1335 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1336 END IF 1337 ELSE 1338 misfdep(ji,jj)= misfdep(ji,jj) - 1 1339 risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat ) 1340 END IF 1341 ENDIF 1342 END DO 1343 END DO 1344 1345 ! point V mbathy(ji,jj) == misfdep(ji,jj+1) 1346 DO jj = 1, jpjm1 1347 DO ji = 1, jpim1 1348 IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1349 IF ( .NOT. ln_iscpl ) THEN 1350 zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) & 1351 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat ))) 1352 zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) & 1353 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat ))) 1354 IF (zbathydiff <= zrisfdepdiff) THEN 1355 mbathy(ji,jj) = mbathy(ji,jj) + 1 1356 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat ) 1357 ELSE 1358 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1359 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1360 END IF 1361 ELSE 1362 misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1 1363 risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat ) 1364 END IF 1365 ENDIF 1366 END DO 1367 END DO 1368 1369 IF( lk_mpp ) THEN 1370 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1371 CALL lbc_lnk( zbathy, 'T', 1. ) 1372 misfdep(:,:) = INT( zbathy(:,:) ) 1373 1374 CALL lbc_lnk( risfdep,'T', 1. ) 1375 CALL lbc_lnk( bathy, 'T', 1. ) 1376 1377 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1378 CALL lbc_lnk( zbathy, 'T', 1. ) 1379 mbathy(:,:) = INT( zbathy(:,:) ) 1380 ENDIF 1381 ! point V misdep(ji,jj) == mbathy(ji,jj+1) 1382 DO jj = 1, jpjm1 1383 DO ji = 1, jpim1 1384 IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN 1385 IF ( .NOT. ln_iscpl ) THEN 1386 zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) & 1387 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ))) 1388 zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) & 1389 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat ))) 1390 IF (zbathydiff <= zrisfdepdiff) THEN 1391 mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1 1392 bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat ) 1393 ELSE 1394 misfdep(ji,jj) = misfdep(ji,jj) - 1 1395 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1396 END IF 1397 ELSE 1398 misfdep(ji,jj) = misfdep(ji,jj) - 1 1399 risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat ) 1400 END IF 1401 ENDIF 1402 END DO 1403 END DO 1404 1405 1406 IF( lk_mpp ) THEN 1407 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1408 CALL lbc_lnk( zbathy, 'T', 1. ) 1409 misfdep(:,:) = INT( zbathy(:,:) ) 1410 1411 CALL lbc_lnk( risfdep,'T', 1. ) 1412 CALL lbc_lnk( bathy, 'T', 1. ) 1413 1414 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1415 CALL lbc_lnk( zbathy, 'T', 1. ) 1416 mbathy(:,:) = INT( zbathy(:,:) ) 1417 ENDIF 1418 1419 ! point U mbathy(ji,jj) == misfdep(ji,jj+1) 1420 DO jj = 1, jpjm1 1421 DO ji = 1, jpim1 1422 IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN 1423 IF ( .NOT. ln_iscpl ) THEN 1424 zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) & 1425 & + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat ))) 1426 zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) & 1427 & - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat ))) 1428 IF (zbathydiff <= zrisfdepdiff) THEN 1429 mbathy(ji,jj) = mbathy(ji,jj) + 1 1430 bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat ) 1431 ELSE 1432 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1433 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1434 END IF 1435 ELSE 1436 misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1 1437 risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat ) 1438 ENDIF 1439 ENDIF 1440 ENDDO 1441 ENDDO 1442 1443 IF( lk_mpp ) THEN 1444 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1445 CALL lbc_lnk( zbathy, 'T', 1. ) 1446 misfdep(:,:) = INT( zbathy(:,:) ) 1447 1448 CALL lbc_lnk( risfdep,'T', 1. ) 1449 CALL lbc_lnk( bathy, 'T', 1. ) 1450 1451 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1452 CALL lbc_lnk( zbathy, 'T', 1. ) 1453 mbathy(:,:) = INT( zbathy(:,:) ) 1454 ENDIF 1455 1456 ! point U misfdep(ji,jj) == bathy(ji,jj+1) 1457 DO jj = 1, jpjm1 1458 DO ji = 1, jpim1 1459 IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN 1460 IF ( .NOT. ln_iscpl ) THEN 1461 zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) & 1462 & + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat ))) 1463 zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) & 1464 & - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat ))) 1465 IF (zbathydiff <= zrisfdepdiff) THEN 1466 mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1 1467 bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat ) 1468 ELSE 1469 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1470 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1471 END IF 1472 ELSE 1473 misfdep(ji,jj) = misfdep(ji ,jj) - 1 1474 risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat ) 1475 ENDIF 1476 ENDIF 1477 ENDDO 1478 ENDDO 1479 1480 IF( lk_mpp ) THEN 1481 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1482 CALL lbc_lnk( zbathy, 'T', 1. ) 1483 misfdep(:,:) = INT( zbathy(:,:) ) 1484 1485 CALL lbc_lnk( risfdep,'T', 1. ) 1486 CALL lbc_lnk( bathy, 'T', 1. ) 1487 1488 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1489 CALL lbc_lnk( zbathy, 'T', 1. ) 1490 mbathy(:,:) = INT( zbathy(:,:) ) 1491 ENDIF 1492 END DO 1493 ! end dig bathy/ice shelf to be compatible 1494 ! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness 1495 DO jl = 1,20 1496 1497 ! remove single point "bay" on isf coast line in the ice shelf draft' 1498 DO jk = 2, jpk 1499 WHERE (misfdep==0) misfdep=jpk 1500 zmask=0._wp 1501 WHERE (misfdep <= jk) zmask=1 1502 DO jj = 2, jpjm1 1503 DO ji = 2, jpim1 1504 IF (misfdep(ji,jj) == jk) THEN 1505 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1506 IF (ibtest <= 1) THEN 1507 risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1 1508 IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk 1509 END IF 1510 END IF 1511 END DO 1512 END DO 1513 END DO 1514 WHERE (misfdep==jpk) 1515 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1516 END WHERE 1517 IF( lk_mpp ) THEN 1518 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1519 CALL lbc_lnk( zbathy, 'T', 1. ) 1520 misfdep(:,:) = INT( zbathy(:,:) ) 1521 1522 CALL lbc_lnk( risfdep,'T', 1. ) 1523 CALL lbc_lnk( bathy, 'T', 1. ) 1524 1525 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1526 CALL lbc_lnk( zbathy, 'T', 1. ) 1527 mbathy(:,:) = INT( zbathy(:,:) ) 1528 ENDIF 1529 1530 ! remove single point "bay" on bathy coast line beneath an ice shelf' 1531 DO jk = jpk,1,-1 1532 zmask=0._wp 1533 WHERE (mbathy >= jk ) zmask=1 1534 DO jj = 2, jpjm1 1535 DO ji = 2, jpim1 1536 IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN 1537 ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1) 1538 IF (ibtest <= 1) THEN 1539 bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1 1540 IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0 1541 END IF 1542 END IF 1543 END DO 1544 END DO 1545 END DO 1546 WHERE (mbathy==0) 1547 misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp 1548 END WHERE 1549 IF( lk_mpp ) THEN 1550 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1551 CALL lbc_lnk( zbathy, 'T', 1. ) 1552 misfdep(:,:) = INT( zbathy(:,:) ) 1553 1554 CALL lbc_lnk( risfdep,'T', 1. ) 1555 CALL lbc_lnk( bathy, 'T', 1. ) 1556 1557 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1558 CALL lbc_lnk( zbathy, 'T', 1. ) 1559 mbathy(:,:) = INT( zbathy(:,:) ) 1560 ENDIF 1561 1562 ! fill hole in ice shelf 1563 zmisfdep = misfdep 1564 zrisfdep = risfdep 1565 WHERE (zmisfdep <= 1._wp) zmisfdep=jpk 1566 DO jj = 2, jpjm1 1567 DO ji = 2, jpim1 1568 ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj ) 1569 ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1) 1570 IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk 1571 IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk 1572 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk 1573 IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk 1574 ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1575 IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN 1576 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp 1577 END IF 1578 IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN 1579 misfdep(ji,jj) = ibtest 1580 risfdep(ji,jj) = gdepw_1d(ibtest) 1581 ENDIF 1582 ENDDO 1583 ENDDO 1584 1585 IF( lk_mpp ) THEN 1586 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1587 CALL lbc_lnk( zbathy, 'T', 1. ) 1588 misfdep(:,:) = INT( zbathy(:,:) ) 1589 1590 CALL lbc_lnk( risfdep, 'T', 1. ) 1591 CALL lbc_lnk( bathy, 'T', 1. ) 1592 1593 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1594 CALL lbc_lnk( zbathy, 'T', 1. ) 1595 mbathy(:,:) = INT( zbathy(:,:) ) 1596 ENDIF 1597 ! 1598 !! fill hole in bathymetry 1599 zmbathy (:,:)=mbathy (:,:) 1600 DO jj = 2, jpjm1 1601 DO ji = 2, jpim1 1602 ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj ) 1603 ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1) 1604 IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0 1605 IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0 1606 IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0 1607 IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0 1608 ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1) 1609 IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN 1610 mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ; 1611 END IF 1612 IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN 1613 mbathy(ji,jj) = ibtest 1614 bathy(ji,jj) = gdepw_1d(ibtest+1) 1615 ENDIF 1616 END DO 1617 END DO 1618 IF( lk_mpp ) THEN 1619 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1620 CALL lbc_lnk( zbathy, 'T', 1. ) 1621 misfdep(:,:) = INT( zbathy(:,:) ) 1622 1623 CALL lbc_lnk( risfdep, 'T', 1. ) 1624 CALL lbc_lnk( bathy, 'T', 1. ) 1625 1626 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1627 CALL lbc_lnk( zbathy, 'T', 1. ) 1628 mbathy(:,:) = INT( zbathy(:,:) ) 1629 ENDIF 1630 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1631 DO jj = 1, jpjm1 1632 DO ji = 1, jpim1 1633 IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1634 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1635 END IF 1636 END DO 1637 END DO 1638 IF( lk_mpp ) THEN 1639 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1640 CALL lbc_lnk( zbathy, 'T', 1. ) 1641 misfdep(:,:) = INT( zbathy(:,:) ) 1642 1643 CALL lbc_lnk( risfdep, 'T', 1. ) 1644 CALL lbc_lnk( bathy, 'T', 1. ) 1645 1646 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1647 CALL lbc_lnk( zbathy, 'T', 1. ) 1648 mbathy(:,:) = INT( zbathy(:,:) ) 1649 ENDIF 1650 ! if not compatible after all check (ie U point water column less than 2 cells), mask U 1651 DO jj = 1, jpjm1 1652 DO ji = 1, jpim1 1653 IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN 1654 mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ; 1655 END IF 1656 END DO 1657 END DO 1658 IF( lk_mpp ) THEN 1659 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1660 CALL lbc_lnk( zbathy, 'T', 1. ) 1661 misfdep(:,:) = INT( zbathy(:,:) ) 1662 1663 CALL lbc_lnk( risfdep,'T', 1. ) 1664 CALL lbc_lnk( bathy, 'T', 1. ) 1665 1666 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1667 CALL lbc_lnk( zbathy, 'T', 1. ) 1668 mbathy(:,:) = INT( zbathy(:,:) ) 1669 ENDIF 1670 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1671 DO jj = 1, jpjm1 1672 DO ji = 1, jpi 1673 IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1674 mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ; 1675 END IF 1676 END DO 1677 END DO 1678 IF( lk_mpp ) THEN 1679 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1680 CALL lbc_lnk( zbathy, 'T', 1. ) 1681 misfdep(:,:) = INT( zbathy(:,:) ) 1682 1683 CALL lbc_lnk( risfdep,'T', 1. ) 1684 CALL lbc_lnk( bathy, 'T', 1. ) 1685 1686 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1687 CALL lbc_lnk( zbathy, 'T', 1. ) 1688 mbathy(:,:) = INT( zbathy(:,:) ) 1689 ENDIF 1690 ! if not compatible after all check (ie V point water column less than 2 cells), mask V 1691 DO jj = 1, jpjm1 1692 DO ji = 1, jpi 1693 IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN 1694 mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ; 1695 END IF 1696 END DO 1697 END DO 1698 IF( lk_mpp ) THEN 1699 zbathy(:,:) = FLOAT( misfdep(:,:) ) 1700 CALL lbc_lnk( zbathy, 'T', 1. ) 1701 misfdep(:,:) = INT( zbathy(:,:) ) 1702 1703 CALL lbc_lnk( risfdep,'T', 1. ) 1704 CALL lbc_lnk( bathy, 'T', 1. ) 1705 1706 zbathy(:,:) = FLOAT( mbathy(:,:) ) 1707 CALL lbc_lnk( zbathy, 'T', 1. ) 1708 mbathy(:,:) = INT( zbathy(:,:) ) 1709 ENDIF 1710 ! if not compatible after all check, mask T 1711 DO jj = 1, jpj 1712 DO ji = 1, jpi 1713 IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN 1714 misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ; 1715 END IF 1716 END DO 1717 END DO 1718 1719 WHERE (mbathy(:,:) == 1) 1720 mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp 1721 END WHERE 1722 END DO 1723 ! end check compatibility ice shelf/bathy 1724 ! remove very shallow ice shelf (less than ~ 10m if 75L) 1725 WHERE (risfdep(:,:) <= 10._wp) 1726 misfdep = 1; risfdep = 0.0_wp; 1727 END WHERE 1728 1729 IF( icompt == 0 ) THEN 1730 IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry' 1731 ELSE 1732 IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry' 1733 ENDIF 1734 1735 ! compute scale factor and depth at T- and W- points 1736 DO jj = 1, jpj 1737 DO ji = 1, jpi 1738 ik = mbathy(ji,jj) 1739 IF( ik > 0 ) THEN ! ocean point only 1740 ! max ocean level case 1741 IF( ik == jpkm1 ) THEN 1742 zdepwp = bathy(ji,jj) 1743 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1744 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1745 e3t_0(ji,jj,ik ) = ze3tp 1746 e3t_0(ji,jj,ik+1) = ze3tp 1747 e3w_0(ji,jj,ik ) = ze3wp 1748 e3w_0(ji,jj,ik+1) = ze3tp 1749 gdepw_0(ji,jj,ik+1) = zdepwp 1750 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1751 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1752 ! 1753 ELSE ! standard case 1754 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1755 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1756 ENDIF 1757 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1758 !gm Bug? check the gdepw_1d 1759 ! ... on ik 1760 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1761 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1762 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1763 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1764 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1765 ! ... on ik+1 1766 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1767 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1768 ENDIF 1769 ENDIF 1770 END DO 1771 END DO 1772 ! 1773 it = 0 1774 DO jj = 1, jpj 1775 DO ji = 1, jpi 1776 ik = mbathy(ji,jj) 1777 IF( ik > 0 ) THEN ! ocean point only 1778 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1779 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1780 ! test 1781 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1782 IF( zdiff <= 0._wp .AND. lwp ) THEN 1783 it = it + 1 1784 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1785 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1786 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1787 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1788 ENDIF 1789 ENDIF 1790 END DO 1791 END DO 1792 ! 1793 ! (ISF) Definition of e3t, u, v, w for ISF case 1794 DO jj = 1, jpj 1795 DO ji = 1, jpi 1796 ik = misfdep(ji,jj) 1797 IF( ik > 1 ) THEN ! ice shelf point only 1798 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1799 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1800 !gm Bug? check the gdepw_0 1801 ! ... on ik 1802 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1803 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1804 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1805 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1806 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1807 1808 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1809 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1810 ENDIF 1811 ! ... on ik / ik-1 1812 e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1813 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1814 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1815 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1816 ENDIF 1817 END DO 1818 END DO 1819 1820 it = 0 1821 DO jj = 1, jpj 1822 DO ji = 1, jpi 1823 ik = misfdep(ji,jj) 1824 IF( ik > 1 ) THEN ! ice shelf point only 1825 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1826 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1827 ! test 1828 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1829 IF( zdiff <= 0. .AND. lwp ) THEN 1830 it = it + 1 1831 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1832 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1833 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1834 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1835 ENDIF 1836 ENDIF 1837 END DO 1838 END DO 1839 1840 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1841 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1842 ! 1843 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1844 ! 1845 END SUBROUTINE zgr_isf 1846 1847 1848 SUBROUTINE zgr_sco 1849 !!---------------------------------------------------------------------- 1850 !! *** ROUTINE zgr_sco *** 1851 !! 1852 !! ** Purpose : define the s-coordinate system 1853 !! 1854 !! ** Method : s-coordinate 1855 !! The depth of model levels is defined as the product of an 1856 !! analytical function by the local bathymetry, while the vertical 1857 !! scale factors are defined as the product of the first derivative 1858 !! of the analytical function by the bathymetry. 1859 !! (this solution save memory as depth and scale factors are not 1860 !! 3d fields) 1861 !! - Read bathymetry (in meters) at t-point and compute the 1862 !! bathymetry at u-, v-, and f-points. 1863 !! hbatu = mi( hbatt ) 1864 !! hbatv = mj( hbatt ) 1865 !! hbatf = mi( mj( hbatt ) ) 1866 !! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical 1867 !! function and its derivative given as function. 1868 !! z_gsigt(k) = fssig (k ) 1869 !! z_gsigw(k) = fssig (k-0.5) 1870 !! z_esigt(k) = fsdsig(k ) 1871 !! z_esigw(k) = fsdsig(k-0.5) 1872 !! Three options for stretching are give, and they can be modified 1873 !! following the users requirements. Nevertheless, the output as 1874 !! well as the way to compute the model levels and scale factors 1875 !! must be respected in order to insure second order accuracy 1876 !! schemes. 1877 !! 1878 !! The three methods for stretching available are: 1879 !! 1880 !! s_sh94 (Song and Haidvogel 1994) 1881 !! a sinh/tanh function that allows sigma and stretched sigma 1882 !! 1883 !! s_sf12 (Siddorn and Furner 2012?) 1884 !! allows the maintenance of fixed surface and or 1885 !! bottom cell resolutions (cf. geopotential coordinates) 1886 !! within an analytically derived stretched S-coordinate framework. 1887 !! 1888 !! s_tanh (Madec et al 1996) 1889 !! a cosh/tanh function that gives stretched coordinates 1890 !! 1891 !!---------------------------------------------------------------------- 1892 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1893 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers 1894 INTEGER :: ios ! Local integer output status for namelist read 1895 REAL(wp) :: zrmax, ztaper ! temporary scalars 1896 REAL(wp) :: zrfact 1897 ! 1898 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1899 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1900 !! 1901 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1902 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1903 !!---------------------------------------------------------------------- 1904 ! 1905 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1906 ! 1907 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1908 ! 1909 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters 1910 READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901) 1911 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp ) 1912 1913 REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters 1914 READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 ) 1915 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp ) 1916 IF(lwm) WRITE ( numond, namzgr_sco ) 1917 1918 IF(lwp) THEN ! control print 1919 WRITE(numout,*) 1920 WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate' 1921 WRITE(numout,*) '~~~~~~~~~~~' 1922 WRITE(numout,*) ' Namelist namzgr_sco' 1923 WRITE(numout,*) ' stretching coeffs ' 1924 WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max 1925 WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min 1926 WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc 1927 WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax 1928 WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94 1929 WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients' 1930 WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta 1931 WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb 1932 WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb 1933 WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12 1934 WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit 1935 WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients' 1936 WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha 1937 WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold 1938 WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs 1939 WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a 1940 WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b 1941 WRITE(numout,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b' 1942 ENDIF 1943 1944 hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate 1945 hifu(:,:) = rn_sbot_min 1946 hifv(:,:) = rn_sbot_min 1947 hiff(:,:) = rn_sbot_min 1948 1949 ! ! set maximum ocean depth 1950 bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) 1951 1952 IF( .NOT.ln_wd ) THEN 1953 DO jj = 1, jpj 1954 DO ji = 1, jpi 1955 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) ) 1956 END DO 1957 END DO 1958 END IF 1959 ! ! ============================= 1960 ! ! Define the envelop bathymetry (hbatt) 1961 ! ! ============================= 1962 ! use r-value to create hybrid coordinates 1963 zenv(:,:) = bathy(:,:) 1964 ! 1965 IF( .NOT.ln_wd ) THEN 1966 ! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing 1967 DO jj = 1, jpj 1968 DO ji = 1, jpi 1969 IF( bathy(ji,jj) == 0._wp ) THEN 1970 iip1 = MIN( ji+1, jpi ) 1971 ijp1 = MIN( jj+1, jpj ) 1972 iim1 = MAX( ji-1, 1 ) 1973 ijm1 = MAX( jj-1, 1 ) 1974 !!gm BUG fix see ticket #1617 1975 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1976 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1977 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) & 1978 & zenv(ji,jj) = rn_sbot_min 1979 !!gm 1980 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1981 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1982 !!gm zenv(ji,jj) = rn_sbot_min 1983 !!gm ENDIF 1984 !!gm end 1985 ENDIF 1986 END DO 1987 END DO 1988 END IF 1989 1990 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 1991 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 1992 ! 1993 ! smooth the bathymetry (if required) 1994 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 1995 scobot(:,:) = bathy(:,:) ! ocean bottom depth 1996 ! 1997 jl = 0 1998 zrmax = 1._wp 1999 ! 2000 ! 2001 ! set scaling factor used in reducing vertical gradients 2002 zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax ) 2003 ! 2004 ! initialise temporary evelope depth arrays 2005 ztmpi1(:,:) = zenv(:,:) 2006 ztmpi2(:,:) = zenv(:,:) 2007 ztmpj1(:,:) = zenv(:,:) 2008 ztmpj2(:,:) = zenv(:,:) 2009 ! 2010 ! initialise temporary r-value arrays 2011 zri(:,:) = 1._wp 2012 zrj(:,:) = 1._wp 2013 ! ! ================ ! 2014 DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop ! 2015 ! ! ================ ! 2016 jl = jl + 1 2017 zrmax = 0._wp 2018 ! we set zrmax from previous r-values (zri and zrj) first 2019 ! if set after current r-value calculation (as previously) 2020 ! we could exit DO WHILE prematurely before checking r-value 2021 ! of current zenv 2022 DO jj = 1, nlcj 2023 DO ji = 1, nlci 2024 zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) ) 2025 END DO 2026 END DO 2027 zri(:,:) = 0._wp 2028 zrj(:,:) = 0._wp 2029 DO jj = 1, nlcj 2030 DO ji = 1, nlci 2031 iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi) 2032 ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj) 2033 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN 2034 zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) ) 2035 END IF 2036 IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN 2037 zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) ) 2038 END IF 2039 IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact 2040 IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact 2041 IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact 2042 IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact 2043 END DO 2044 END DO 2045 IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain 2046 ! 2047 IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax 2048 ! 2049 DO jj = 1, nlcj 2050 DO ji = 1, nlci 2051 zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) ) 2052 END DO 2053 END DO 2054 ! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero 2055 CALL lbc_lnk( zenv, 'T', 1._wp, 'no0' ) 2056 ! ! ================ ! 2057 END DO ! End loop ! 2058 ! ! ================ ! 2059 DO jj = 1, jpj 2060 DO ji = 1, jpi 2061 zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings 2062 END DO 2063 END DO 2064 ! 2065 ! Envelope bathymetry saved in hbatt 2066 hbatt(:,:) = zenv(:,:) 2067 IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN 2068 CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' ) 2069 DO jj = 1, jpj 2070 DO ji = 1, jpi 2071 ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp ) 2072 hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper ) 2073 END DO 2074 END DO 2075 ENDIF 2076 ! 2077 ! ! ============================== 2078 ! ! hbatu, hbatv, hbatf fields 2079 ! ! ============================== 2080 IF(lwp) THEN 2081 WRITE(numout,*) 2082 IF( .NOT.ln_wd ) THEN 2083 WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min 2084 ELSE 2085 WRITE(numout,*) ' zgr_sco: minimum positive depth of the envelop topography set to : ', rn_sbot_min 2086 WRITE(numout,*) ' zgr_sco: minimum negative depth of the envelop topography set to : ', -rn_wdld 2087 ENDIF 2088 ENDIF 2089 hbatu(:,:) = rn_sbot_min 2090 hbatv(:,:) = rn_sbot_min 2091 hbatf(:,:) = rn_sbot_min 2092 DO jj = 1, jpjm1 2093 DO ji = 1, jpim1 ! NO vector opt. 2094 hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) ) 2095 hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) ) 2096 hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) & 2097 & + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) ) 2098 END DO 2099 END DO 2100 2101 IF( ln_wd ) THEN !avoid the zero depth on T- (U-,V-,F-) points 2102 DO jj = 1, jpj 2103 DO ji = 1, jpi 2104 IF(ABS(hbatt(ji,jj)) < rn_wdmin1) & 2105 & hbatt(ji,jj) = SIGN(1._wp, hbatt(ji,jj)) * rn_wdmin1 2106 IF(ABS(hbatu(ji,jj)) < rn_wdmin1) & 2107 & hbatu(ji,jj) = SIGN(1._wp, hbatu(ji,jj)) * rn_wdmin1 2108 IF(ABS(hbatv(ji,jj)) < rn_wdmin1) & 2109 & hbatv(ji,jj) = SIGN(1._wp, hbatv(ji,jj)) * rn_wdmin1 2110 IF(ABS(hbatf(ji,jj)) < rn_wdmin1) & 2111 & hbatf(ji,jj) = SIGN(1._wp, hbatf(ji,jj)) * rn_wdmin1 2112 END DO 2113 END DO 2114 END IF 2115 ! 2116 ! Apply lateral boundary condition 2117 !!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL 2118 zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp ) 2119 DO jj = 1, jpj 2120 DO ji = 1, jpi 2121 IF( hbatu(ji,jj) == 0._wp ) THEN 2122 !No worries about the following line when ln_wd == .true. 2123 IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min 2124 IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj) 2125 ENDIF 2126 END DO 2127 END DO 2128 zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp ) 2129 DO jj = 1, jpj 2130 DO ji = 1, jpi 2131 IF( hbatv(ji,jj) == 0._wp ) THEN 2132 IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min 2133 IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj) 2134 ENDIF 2135 END DO 2136 END DO 2137 zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp ) 2138 DO jj = 1, jpj 2139 DO ji = 1, jpi 2140 IF( hbatf(ji,jj) == 0._wp ) THEN 2141 IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min 2142 IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj) 2143 ENDIF 2144 END DO 2145 END DO 2146 2147 !!bug: key_helsinki a verifer 2148 IF( .NOT.ln_wd ) THEN 2149 hift(:,:) = MIN( hift(:,:), hbatt(:,:) ) 2150 hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) ) 2151 hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) ) 2152 hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) ) 2153 END IF 2154 2155 IF( nprint == 1 .AND. lwp ) THEN 2156 WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), & 2157 & ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) ) 2158 WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), & 2159 & ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) ) 2160 WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), & 2161 & ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) ) 2162 WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), & 2163 & ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) ) 2164 ENDIF 2165 !! helsinki 2166 2167 ! ! ======================= 2168 ! ! s-ccordinate fields (gdep., e3.) 2169 ! ! ======================= 2170 ! 2171 ! non-dimensional "sigma" for model level depth at w- and t-levels 2172 2173 2174 !======================================================================== 2175 ! Song and Haidvogel 1994 (ln_s_sh94=T) 2176 ! Siddorn and Furner 2012 (ln_sf12=T) 2177 ! or tanh function (both false) 2178 !======================================================================== 2179 IF ( ln_s_sh94 ) THEN 2180 CALL s_sh94() 2181 ELSE IF ( ln_s_sf12 ) THEN 2182 CALL s_sf12() 2183 ELSE 2184 CALL s_tanh() 2185 ENDIF 2186 2187 CALL lbc_lnk( e3t_0 , 'T', 1._wp ) 2188 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) 2189 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) 2190 CALL lbc_lnk( e3f_0 , 'F', 1._wp ) 2191 CALL lbc_lnk( e3w_0 , 'W', 1._wp ) 2192 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2193 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2194 ! 2195 IF( .NOT.ln_wd ) THEN 2196 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2197 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2198 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2199 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2200 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2201 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2202 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2203 END IF 2204 2205 #if defined key_agrif 2206 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2207 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2208 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2209 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2210 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2211 ENDIF 2212 #endif 2213 2214 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2215 !!gm and only that !!!!! 2216 !!gm THIS should be removed from here ! 2217 gdept_n(:,:,:) = gdept_0(:,:,:) 2218 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2219 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2220 e3t_n (:,:,:) = e3t_0 (:,:,:) 2221 e3u_n (:,:,:) = e3u_0 (:,:,:) 2222 e3v_n (:,:,:) = e3v_0 (:,:,:) 2223 e3f_n (:,:,:) = e3f_0 (:,:,:) 2224 e3w_n (:,:,:) = e3w_0 (:,:,:) 2225 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2226 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2227 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2228 !! gm end 2229 !! 2230 ! HYBRID : 2231 DO jj = 1, jpj 2232 DO ji = 1, jpi 2233 DO jk = 1, jpkm1 2234 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2235 END DO 2236 IF( ln_wd ) THEN 2237 IF( scobot(ji,jj) <= -(rn_wdld - rn_wdmin2)) THEN 2238 mbathy(ji,jj) = 0 2239 ELSEIF(scobot(ji,jj) <= rn_wdmin1) THEN 2240 mbathy(ji,jj) = 2 2241 ENDIF 2242 ELSE 2243 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2244 ENDIF 2245 END DO 2246 END DO 2247 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 2248 & ' MAX ', MAXVAL( mbathy(:,:) ) 2249 2250 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 2251 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2252 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2253 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2254 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2255 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2256 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2257 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2258 2259 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2260 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2261 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2262 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2263 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2264 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2265 ENDIF 2266 ! END DO 2267 IF(lwp) THEN ! selected vertical profiles 2268 WRITE(numout,*) 2269 WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1) 2270 WRITE(numout,*) ' ~~~~~~ --------------------' 2271 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 2272 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), & 2273 & e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk ) 2274 DO jj = mj0(20), mj1(20) 2275 DO ji = mi0(20), mi1(20) 2276 WRITE(numout,*) 2277 WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 2278 WRITE(numout,*) ' ~~~~~~ --------------------' 2279 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 2280 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 2281 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 2282 END DO 2283 END DO 2284 DO jj = mj0(74), mj1(74) 2285 DO ji = mi0(100), mi1(100) 2286 WRITE(numout,*) 2287 WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj) 2288 WRITE(numout,*) ' ~~~~~~ --------------------' 2289 WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')") 2290 WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), & 2291 & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk ) 2292 END DO 2293 END DO 2294 ENDIF 2295 ! 2296 !================================================================================ 2297 ! check the coordinate makes sense 2298 !================================================================================ 2299 DO ji = 1, jpi 2300 DO jj = 1, jpj 2301 ! 2302 IF( hbatt(ji,jj) > 0._wp) THEN 2303 DO jk = 1, mbathy(ji,jj) 2304 ! check coordinate is monotonically increasing 2305 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2306 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2307 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2308 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2309 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2310 CALL ctl_stop( ctmp1 ) 2311 ENDIF 2312 ! and check it has never gone negative 2313 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2314 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2315 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2316 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2317 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2318 CALL ctl_stop( ctmp1 ) 2319 ENDIF 2320 ! and check it never exceeds the total depth 2321 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2322 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2323 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2324 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2325 CALL ctl_stop( ctmp1 ) 2326 ENDIF 2327 END DO 2328 ! 2329 DO jk = 1, mbathy(ji,jj)-1 2330 ! and check it never exceeds the total depth 2331 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2332 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2333 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2334 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2335 CALL ctl_stop( ctmp1 ) 2336 ENDIF 2337 END DO 2338 ENDIF 2339 END DO 2340 END DO 2341 ! 2342 CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 2343 ! 2344 IF( nn_timing == 1 ) CALL timing_stop('zgr_sco') 2345 ! 2346 END SUBROUTINE zgr_sco 2347 2348 2349 SUBROUTINE s_sh94() 2350 !!---------------------------------------------------------------------- 2351 !! *** ROUTINE s_sh94 *** 2352 !! 2353 !! ** Purpose : stretch the s-coordinate system 2354 !! 2355 !! ** Method : s-coordinate stretch using the Song and Haidvogel 1994 2356 !! mixed S/sigma coordinate 2357 !! 2358 !! Reference : Song and Haidvogel 1994. 2359 !!---------------------------------------------------------------------- 2360 INTEGER :: ji, jj, jk ! dummy loop argument 2361 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2362 REAL(wp) :: ztmpu, ztmpv, ztmpf 2363 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2364 ! 2365 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2366 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2367 !!---------------------------------------------------------------------- 2368 2369 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2370 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2371 2372 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 2373 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 2374 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2375 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2376 ! 2377 DO ji = 1, jpi 2378 DO jj = 1, jpj 2379 ! 2380 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2381 DO jk = 1, jpk 2382 z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb ) 2383 z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb ) 2384 END DO 2385 ELSE ! shallow water, uniform sigma 2386 DO jk = 1, jpk 2387 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp) 2388 z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp) 2389 END DO 2390 ENDIF 2391 ! 2392 DO jk = 1, jpkm1 2393 z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 2394 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 2395 END DO 2396 z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) ) 2397 z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) ) 2398 ! 2399 ! Coefficients for vertical depth as the sum of e3w scale factors 2400 z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1) 2401 DO jk = 2, jpk 2402 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 2403 END DO 2404 ! 2405 DO jk = 1, jpk 2406 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2407 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2408 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2409 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2410 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2411 END DO 2412 ! 2413 END DO ! for all jj's 2414 END DO ! for all ji's 2415 2416 DO ji = 1, jpim1 2417 DO jj = 1, jpjm1 2418 ! extended for Wetting/Drying case 2419 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2420 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2421 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2422 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2423 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2424 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2425 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2426 DO jk = 1, jpk 2427 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2428 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2429 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2430 ELSE 2431 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2432 & / ztmpu 2433 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2434 & / ztmpu 2435 END IF 2436 2437 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2438 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2439 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2440 ELSE 2441 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2442 & / ztmpv 2443 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2444 & / ztmpv 2445 END IF 2446 2447 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2448 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj ,jk) + z_esigt3(ji+1,jj ,jk) & 2449 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2450 ELSE 2451 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2452 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2453 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2454 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2455 END IF 2456 2457 ! 2458 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2459 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2460 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2461 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2462 ! 2463 e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2464 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2465 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) ) 2466 END DO 2467 END DO 2468 END DO 2469 ! 2470 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2471 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2472 ! 2473 END SUBROUTINE s_sh94 2474 2475 2476 SUBROUTINE s_sf12 2477 !!---------------------------------------------------------------------- 2478 !! *** ROUTINE s_sf12 *** 2479 !! 2480 !! ** Purpose : stretch the s-coordinate system 2481 !! 2482 !! ** Method : s-coordinate stretch using the Siddorn and Furner 2012? 2483 !! mixed S/sigma/Z coordinate 2484 !! 2485 !! This method allows the maintenance of fixed surface and or 2486 !! bottom cell resolutions (cf. geopotential coordinates) 2487 !! within an analytically derived stretched S-coordinate framework. 2488 !! 2489 !! 2490 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2491 !!---------------------------------------------------------------------- 2492 INTEGER :: ji, jj, jk ! dummy loop argument 2493 REAL(wp) :: zsmth ! smoothing around critical depth 2494 REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space 2495 REAL(wp) :: ztmpu, ztmpv, ztmpf 2496 REAL(wp) :: ztmpu1, ztmpv1, ztmpf1 2497 ! 2498 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2499 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2500 !!---------------------------------------------------------------------- 2501 ! 2502 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2503 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2504 2505 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp 2506 z_esigt3 = 0._wp ; z_esigw3 = 0._wp 2507 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2508 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2509 2510 DO ji = 1, jpi 2511 DO jj = 1, jpj 2512 2513 IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma 2514 2515 zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,. 2516 ! could be changed by users but care must be taken to do so carefully 2517 zzb = 1.0_wp-(zzb/hbatt(ji,jj)) 2518 2519 zzs = rn_zs / hbatt(ji,jj) 2520 2521 IF (rn_efold /= 0.0_wp) THEN 2522 zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold ) 2523 ELSE 2524 zsmth = 1.0_wp 2525 ENDIF 2526 2527 DO jk = 1, jpk 2528 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 2529 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp) 2530 ENDDO 2531 z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth ) 2532 z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth ) 2533 2534 ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma 2535 2536 DO jk = 1, jpk 2537 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp) 2538 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp) 2539 END DO 2540 2541 ELSE ! shallow water, z coordinates 2542 2543 DO jk = 1, jpk 2544 z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 2545 z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj)) 2546 END DO 2547 2548 ENDIF 2549 2550 DO jk = 1, jpkm1 2551 z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk) 2552 z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk) 2553 END DO 2554 z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 )) 2555 z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk)) 2556 2557 ! Coefficients for vertical depth as the sum of e3w scale factors 2558 z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1) 2559 DO jk = 2, jpk 2560 z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk) 2561 END DO 2562 2563 DO jk = 1, jpk 2564 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2565 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2566 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2567 END DO 2568 2569 ENDDO ! for all jj's 2570 ENDDO ! for all ji's 2571 2572 DO ji=1,jpi-1 2573 DO jj=1,jpj-1 2574 2575 ! extend to suit for Wetting/Drying case 2576 ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj) 2577 ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1) 2578 ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) 2579 ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj) 2580 ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1) 2581 ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * & 2582 & MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) 2583 DO jk = 1, jpk 2584 IF( ln_wd .AND. (ztmpu1 < 0._wp.OR.ABS(ztmpu) < rn_wdmin1) ) THEN 2585 z_esigtu3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) ) 2586 z_esigwu3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji+1,jj,jk) ) 2587 ELSE 2588 z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2589 & / ztmpu 2590 z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2591 & / ztmpu 2592 END IF 2593 2594 IF( ln_wd .AND. (ztmpv1 < 0._wp.OR.ABS(ztmpv) < rn_wdmin1) ) THEN 2595 z_esigtv3(ji,jj,jk) = 0.5_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji,jj+1,jk) ) 2596 z_esigwv3(ji,jj,jk) = 0.5_wp * ( z_esigw3(ji,jj,jk) + z_esigw3(ji,jj+1,jk) ) 2597 ELSE 2598 z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2599 & / ztmpv 2600 z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2601 & / ztmpv 2602 END IF 2603 2604 IF( ln_wd .AND. (ztmpf1 < 0._wp.OR.ABS(ztmpf) < rn_wdmin1) ) THEN 2605 z_esigtf3(ji,jj,jk) = 0.25_wp * ( z_esigt3(ji,jj,jk) + z_esigt3(ji+1,jj,jk) & 2606 & + z_esigt3(ji,jj+1,jk) + z_esigt3(ji+1,jj+1,jk) ) 2607 ELSE 2608 z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2609 & + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2610 & + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2611 & + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf 2612 END IF 2613 2614 ! Code prior to wetting and drying option (for reference) 2615 !z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) & 2616 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2617 ! 2618 !z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) & 2619 ! /( hbatt(ji,jj)+hbatt(ji+1,jj) ) 2620 ! 2621 !z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) & 2622 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2623 ! 2624 !z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) & 2625 ! /( hbatt(ji,jj)+hbatt(ji,jj+1) ) 2626 ! 2627 !z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) & 2628 ! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) & 2629 ! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) & 2630 ! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) & 2631 ! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) ) 2632 2633 e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk) 2634 e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk) 2635 e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk) 2636 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2637 ! 2638 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2639 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2640 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) 2641 END DO 2642 2643 ENDDO 2644 ENDDO 2645 ! 2646 CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.) 2647 CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.) 2648 CALL lbc_lnk(e3w_0 ,'T',1.) 2649 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2650 ! 2651 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2652 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2653 ! 2654 END SUBROUTINE s_sf12 2655 2656 2657 SUBROUTINE s_tanh() 2658 !!---------------------------------------------------------------------- 2659 !! *** ROUTINE s_tanh*** 2660 !! 2661 !! ** Purpose : stretch the s-coordinate system 2662 !! 2663 !! ** Method : s-coordinate stretch 2664 !! 2665 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2666 !!---------------------------------------------------------------------- 2667 INTEGER :: ji, jj, jk ! dummy loop argument 2668 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2669 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2670 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2671 !!---------------------------------------------------------------------- 2672 2673 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2674 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2675 2676 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp 2677 z_esigt = 0._wp ; z_esigw = 0._wp 2678 2679 DO jk = 1, jpk 2680 z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp ) 2681 z_gsigt(jk) = -fssig( REAL(jk,wp) ) 2682 END DO 2683 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk) 2684 ! 2685 ! Coefficients for vertical scale factors at w-, t- levels 2686 !!gm bug : define it from analytical function, not like juste bellow.... 2687 !!gm or betteroffer the 2 possibilities.... 2688 DO jk = 1, jpkm1 2689 z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk) 2690 z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk) 2691 END DO 2692 z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) ) 2693 z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) ) 2694 ! 2695 ! Coefficients for vertical depth as the sum of e3w scale factors 2696 z_gsi3w(1) = 0.5_wp * z_esigw(1) 2697 DO jk = 2, jpk 2698 z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk) 2699 END DO 2700 !!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays) 2701 DO jk = 1, jpk 2702 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2703 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2704 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2705 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2706 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2707 END DO 2708 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) 2709 DO jj = 1, jpj 2710 DO ji = 1, jpi 2711 DO jk = 1, jpk 2712 e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 2713 e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 2714 e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 2715 e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) ) 2716 ! 2717 e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) ) 2718 e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) ) 2719 e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) ) 2720 END DO 2721 END DO 2722 END DO 2723 ! 2724 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2725 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2726 ! 2727 END SUBROUTINE s_tanh 2728 2729 2730 FUNCTION fssig( pk ) RESULT( pf ) 2731 !!---------------------------------------------------------------------- 2732 !! *** ROUTINE fssig *** 2733 !! 2734 !! ** Purpose : provide the analytical function in s-coordinate 2735 !! 2736 !! ** Method : the function provide the non-dimensional position of 2737 !! T and W (i.e. between 0 and 1) 2738 !! T-points at integer values (between 1 and jpk) 2739 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 2740 !!---------------------------------------------------------------------- 2741 REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate 2742 REAL(wp) :: pf ! sigma value 2743 !!---------------------------------------------------------------------- 2744 ! 2745 pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) & 2746 & - TANH( rn_thetb * rn_theta ) ) & 2747 & * ( COSH( rn_theta ) & 2748 & + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) & 2749 & / ( 2._wp * SINH( rn_theta ) ) 2750 ! 2751 END FUNCTION fssig 2752 2753 2754 FUNCTION fssig1( pk1, pbb ) RESULT( pf1 ) 2755 !!---------------------------------------------------------------------- 2756 !! *** ROUTINE fssig1 *** 2757 !! 2758 !! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate 2759 !! 2760 !! ** Method : the function provides the non-dimensional position of 2761 !! T and W (i.e. between 0 and 1) 2762 !! T-points at integer values (between 1 and jpk) 2763 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 2764 !!---------------------------------------------------------------------- 2765 REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate 2766 REAL(wp), INTENT(in) :: pbb ! Stretching coefficient 2767 REAL(wp) :: pf1 ! sigma value 2768 !!---------------------------------------------------------------------- 2769 ! 2770 IF ( rn_theta == 0 ) then ! uniform sigma 2771 pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 ) 2772 ELSE ! stretched sigma 2773 pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) & 2774 & + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) & 2775 & / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) ) 2776 ENDIF 2777 ! 2778 END FUNCTION fssig1 2779 2780 2781 FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma ) 2782 !!---------------------------------------------------------------------- 2783 !! *** ROUTINE fgamma *** 2784 !! 2785 !! ** Purpose : provide analytical function for the s-coordinate 2786 !! 2787 !! ** Method : the function provides the non-dimensional position of 2788 !! T and W (i.e. between 0 and 1) 2789 !! T-points at integer values (between 1 and jpk) 2790 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 2791 !! 2792 !! This method allows the maintenance of fixed surface and or 2793 !! bottom cell resolutions (cf. geopotential coordinates) 2794 !! within an analytically derived stretched S-coordinate framework. 2795 !! 2796 !! Reference : Siddorn and Furner, in prep 2797 !!---------------------------------------------------------------------- 2798 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2799 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2800 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2801 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2802 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2803 ! 2804 INTEGER :: jk ! dummy loop index 2805 REAL(wp) :: za1,za2,za3 ! local scalar 2806 REAL(wp) :: zn1,zn2 ! - - 2807 REAL(wp) :: za,zb,zx ! - - 2808 !!---------------------------------------------------------------------- 2809 ! 2810 zn1 = 1._wp / REAL( jpkm1, wp ) 2811 zn2 = 1._wp - zn1 2812 ! 2813 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2814 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2815 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2816 ! 2817 za = pzb - za3*(pzs-za1)-za2 2818 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2819 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2820 zx = 1.0_wp-za/2.0_wp-zb 2821 ! 2822 DO jk = 1, jpk 2823 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2824 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2825 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2826 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2827 END DO 2828 ! 2829 END FUNCTION fgamma 314 zk(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk( zk, 'U', 1. ) ; miku(:,:) = MAX( INT( zk(:,:) ), 1 ) 315 zk(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk( zk, 'V', 1. ) ; mikv(:,:) = MAX( INT( zk(:,:) ), 1 ) 316 zk(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk( zk, 'F', 1. ) ; mikf(:,:) = MAX( INT( zk(:,:) ), 1 ) 317 ! 318 zk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk( zk, 'U', 1. ) ; mbku(:,:) = MAX( INT( zk(:,:) ), 1 ) 319 zk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk( zk, 'V', 1. ) ; mbkv(:,:) = MAX( INT( zk(:,:) ), 1 ) 320 ! 321 CALL wrk_dealloc( jpi,jpj, zk ) 322 ! 323 IF( nn_timing == 1 ) CALL timing_stop('zgr_top_bot') 324 ! 325 END SUBROUTINE zgr_top_bot 2830 326 2831 327 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
