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