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