- Timestamp:
- 2018-07-26T09:50:51+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90
r9939 r10001 4 4 !! Ocean : 5 5 !!====================================================================== 6 !! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code 7 !! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 6 !! History : 5.0 ! 2018-07 (G. Madec) Flux Form with Kinetic Energy conservation 7 !! ==>>> here z* and s* only (no z-tilde) 8 9 ! 1- remove z-tilde ==>>> pure z-star (or s-star) 10 ! 2- remove dom_vvl_interpol 11 10 12 !!---------------------------------------------------------------------- 11 13 … … 14 16 !! dom_vvl_sf_nxt : Compute next vertical scale factors 15 17 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid 16 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another17 18 !! dom_vvl_rst : read/write restart file 18 19 !! dom_vvl_ctl : Check the vvl options … … 38 39 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 40 PUBLIC dom_vvl_sf_swp ! called by step.F90 40 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 PUBLIC ssh2e3_before ! ... 42 PUBLIC ssh2e3_now ! ... 41 43 42 44 ! !!* Namelist nam_vvl … … 61 63 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 62 64 65 66 !!gm add 67 !!gm 68 69 70 63 71 !! * Substitutions 64 72 # include "vectopt_loop_substitute.h90" … … 74 82 !! *** FUNCTION dom_vvl_alloc *** 75 83 !!---------------------------------------------------------------------- 84 ! 76 85 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 77 86 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN … … 115 124 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 116 125 !!---------------------------------------------------------------------- 117 INTEGER :: ji, jj, jk118 INTEGER :: ii0, ii1, ij0, ij1119 REAL(wp):: zcoef, z1_Dt120 126 !!---------------------------------------------------------------------- 121 127 ! … … 129 135 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 130 136 ! 131 ! ! Read or initialize e3t_(b/n), te3t_(b/n) and hdiv_lf137 ! ! Read or initialize e3t_(b/n), ssh(b/n) 132 138 CALL dom_vvl_rst( nit000, 'READ' ) 133 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 134 ! 135 ! !== Set of all other vertical scale factors ==! (now and before) 136 ! ! Horizontal interpolation of e3t 137 CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' ) ! from T to U 138 CALL dom_vvl_interpol( e3t_n(:,:,:), e3u_n(:,:,:), 'U' ) 139 CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' ) ! from T to V 140 CALL dom_vvl_interpol( e3t_n(:,:,:), e3v_n(:,:,:), 'V' ) 141 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) ! from U to F 142 ! ! Vertical interpolation of e3t,u,v 143 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n (:,:,:), 'W' ) ! from T to W 144 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b (:,:,:), 'W' ) 145 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) ! from U to UW 146 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 147 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) ! from V to UW 148 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 149 150 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 151 e3t_a(:,:,:) = e3t_n(:,:,:) 152 e3u_a(:,:,:) = e3u_n(:,:,:) 153 e3v_a(:,:,:) = e3v_n(:,:,:) 154 ! 155 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 156 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 157 gdepw_n(:,:,1) = 0.0_wp 158 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 159 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 160 gdepw_b(:,:,1) = 0.0_wp 161 DO jk = 2, jpk ! vertical sum 162 DO jj = 1,jpj 163 DO ji = 1,jpi 164 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 165 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 166 ! ! 0.5 where jk = mikt 167 !!gm ??????? BUG ? gdept_n as well as gde3w_n does not include the thickness of ISF ?? 168 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 169 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 170 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk)) & 171 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk)) 172 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 173 gdepw_b(ji,jj,jk) = gdepw_b(ji,jj,jk-1) + e3t_b(ji,jj,jk-1) 174 gdept_b(ji,jj,jk) = zcoef * ( gdepw_b(ji,jj,jk ) + 0.5 * e3w_b(ji,jj,jk)) & 175 & + (1-zcoef) * ( gdept_b(ji,jj,jk-1) + e3w_b(ji,jj,jk)) 176 END DO 177 END DO 178 END DO 179 ! 180 ! !== thickness of the water column !! (ocean portion only) 181 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 182 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 183 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 184 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 185 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 186 DO jk = 2, jpkm1 187 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 188 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 189 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 190 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 191 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 192 END DO 193 ! 194 ! !== inverse of water column thickness ==! (u- and v- points) 195 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 196 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 197 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 198 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 199 200 ! !== z_tilde coordinate case ==! (Restoring frequencies) 201 IF( ln_vvl_ztilde ) THEN 202 !!gm : idea: add here a READ in a file of custumized restoring frequency 203 ! ! Values in days provided via the namelist 204 ! ! use rsmall to avoid possible division by zero errors with faulty settings 205 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 206 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 207 ! 208 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 209 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 210 frq_rst_hdv(:,:) = 1._wp / rn_Dt 211 ENDIF 212 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 213 z1_Dt = 1._wp / rn_Dt 214 DO jj = 1, jpj 215 DO ji = 1, jpi 216 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 217 IF( ABS(gphit(ji,jj)) >= 6.) THEN 218 ! values outside the equatorial band and transition zone (ztilde) 219 frq_rst_e3t(ji,jj) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400._wp ) 220 frq_rst_hdv(ji,jj) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400._wp ) 221 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 222 ! values inside the equatorial band (ztilde as zstar) 223 frq_rst_e3t(ji,jj) = 0._wp 224 frq_rst_hdv(ji,jj) = z1_Dt 225 ELSE ! transition band (2.5 to 6 degrees N/S) 226 ! ! (linearly transition from z-tilde to z-star) 227 frq_rst_e3t(ji,jj) = 0._wp + ( frq_rst_e3t(ji,jj) - 0._wp ) * 0.5_wp & 228 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp ) ) 229 frq_rst_hdv(ji,jj) = z1_Dt + ( frq_rst_hdv(ji,jj) - z1_Dt ) * 0.5_wp & 230 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) * 180._wp / 3.5_wp ) ) 231 ENDIF 232 END DO 233 END DO 234 IF( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 235 ii0 = 103 ; ii1 = 111 236 ij0 = 128 ; ij1 = 135 ; 237 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0._wp 238 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = z1_Dt 239 ENDIF 240 ENDIF 241 ENDIF 242 ! 243 IF(lwxios) THEN 244 ! define variables in restart file when writing with XIOS 139 ! 140 ! !== Set of all other vertical mesh fields ==! (now and before) 141 ! 142 ! !* BEFORE fields : 143 CALL ssh2e3_before ! set: hu , hv , r1_hu, r1_hv 144 ! ! e3t, e3w, e3u, e3uw, e3v, e3vw 145 ! 146 ! ! set one for all last level to the e3._0 value 147 e3t_b(:,:,jpk) = e3t_0(:,:,jpk) ; e3u_b(:,:,jpk) = e3w_0(:,:,jpk) ; e3v_b(:,:,jpk) = e3v_0(:,:,jpk) 148 e3w_b(:,:,jpk) = e3w_0(:,:,jpk) ; e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk) ; e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk) 149 ! 150 ! !* NOW fields : 151 CALL ssh2e3_now ! set: ht , hu , hv , r1_hu, r1_hv 152 ! ! e3t, e3w, e3u, e3uw, e3v, e3vw, e3f 153 ! ! gdept_n, gdepw_n, gde3w_n 154 ! 155 ! ! set one for all last level to the e3._0 value 156 e3t_n(:,:,jpk) = e3t_0(:,:,jpk) ; e3u_n(:,:,jpk) = e3w_0(:,:,jpk) ; e3v_n(:,:,jpk) = e3v_0(:,:,jpk) 157 e3w_n(:,:,jpk) = e3w_0(:,:,jpk) ; e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk) ; e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk) 158 e3f_n(:,:,jpk) = e3f_0(:,:,jpk) 159 ! 160 ! !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation) 161 e3t_a(:,:,:) = e3t_n(:,:,:) ; e3u_a(:,:,:) = e3u_n(:,:,:) ; e3v_a(:,:,:) = e3v_n(:,:,:) 162 163 !!gm 164 !!gm ===>>>> below: some issues to think about !!! 165 !!gm 166 !!gm fmask definition checked (0 or 1 nothing else) 167 !! in z-tilde or ALE e3._0 should be the time varying fields step forward with an euler scheme 168 !!gm e3w_b & gdept_b are not used except its update in agrif 169 !! and used to compute before slope of surface in dynldf_iso ==>>> remove it !!!! 170 !! NB: in triads on TRA, gdept_n is used !!!! BUG? 171 !!gm e3f_n almost not used ===>>>> verify whether it can be removed or not... 172 !! verify the use of wumask & wvmask mau be replaced by a multiplication by umask(k)*umask(k+1) ??? 173 !! 174 !!gm ISF case to be checked by Pierre Mathiot 175 !! 176 !!gm setting of e3._a for agrif.... TO BE CHECKED 177 178 ! 179 IF(lwxios) THEN ! define variables in restart file when writing with XIOS 245 180 CALL iom_set_rstw_var_active('e3t_b') 246 181 CALL iom_set_rstw_var_active('e3t_n') 247 ! ! ----------------------- !248 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases !249 ! ! ----------------------- !250 CALL iom_set_rstw_var_active('tilde_e3t_b')251 CALL iom_set_rstw_var_active('tilde_e3t_n')252 END IF253 ! ! -------------!254 IF( ln_vvl_ztilde ) THEN ! z_tilde case !255 ! ! ------------ !256 CALL iom_set_rstw_var_active('hdiv_lf')257 ENDIF258 !259 182 ENDIF 260 183 ! … … 287 210 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 288 211 ! 289 INTEGER :: ji, jj, jk ! dummy loop indices 290 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 291 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 292 LOGICAL :: ll_do_bclinic ! local logical 293 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv 294 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t 212 INTEGER :: ji, jj, jk ! dummy loop indices 213 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h 295 214 !!---------------------------------------------------------------------- 296 215 ! … … 305 224 ENDIF 306 225 307 ll_do_bclinic = .TRUE. 308 IF( PRESENT(kcall) ) THEN 309 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 310 ENDIF 311 312 ! ******************************* ! 313 ! After acale factors at t-points ! 314 ! ******************************* ! 315 ! ! --------------------------------------------- ! 316 ! ! z_star coordinate and barotropic z-tilde part ! 317 ! ! --------------------------------------------- ! 318 ! 319 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 320 DO jk = 1, jpkm1 321 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 322 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 323 END DO 324 ! 325 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 326 ! ! ------baroclinic part------ ! 327 ! I - initialization 328 ! ================== 329 330 ! 1 - barotropic divergence 331 ! ------------------------- 332 zhdiv(:,:) = 0._wp 333 zht(:,:) = 0._wp 334 DO jk = 1, jpkm1 335 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 336 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 337 END DO 338 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 339 340 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 341 ! -------------------------------------------------- 342 IF( ln_vvl_ztilde ) THEN 343 IF( kt > nit000 ) THEN 344 DO jk = 1, jpkm1 345 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:) & 346 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 347 END DO 348 ENDIF 349 ENDIF 350 351 ! II - after z_tilde increments of vertical scale factors 352 ! ======================================================= 353 te3t_a(:,:,:) = 0._wp ! te3t_a used to store tendency terms 354 355 ! 1 - High frequency divergence term 356 ! ---------------------------------- 357 IF( ln_vvl_ztilde ) THEN ! z_tilde case 358 DO jk = 1, jpkm1 359 te3t_a(:,:,jk) = te3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 360 END DO 361 ELSE ! layer case 362 DO jk = 1, jpkm1 363 te3t_a(:,:,jk) = te3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 364 END DO 365 ENDIF 366 367 ! 2 - Restoring term (z-tilde case only) 368 ! ------------------ 369 IF( ln_vvl_ztilde ) THEN 370 DO jk = 1, jpk 371 te3t_a(:,:,jk) = te3t_a(:,:,jk) - frq_rst_e3t(:,:) * te3t_b(:,:,jk) 372 END DO 373 ENDIF 374 375 ! 3 - Thickness diffusion term 376 ! ---------------------------- 377 zwu(:,:) = 0._wp 378 zwv(:,:) = 0._wp 379 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 380 DO jj = 1, jpjm1 381 DO ji = 1, fs_jpim1 ! vector opt. 382 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 383 & * ( te3t_b(ji,jj,jk) - te3t_b(ji+1,jj ,jk) ) 384 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 385 & * ( te3t_b(ji,jj,jk) - te3t_b(ji ,jj+1,jk) ) 386 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 387 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 388 END DO 389 END DO 390 END DO 391 DO jj = 1, jpj ! b - correction for last oceanic u-v points 392 DO ji = 1, jpi 393 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 394 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 395 END DO 396 END DO 397 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 398 DO jj = 2, jpjm1 399 DO ji = fs_2, fs_jpim1 ! vector opt. 400 te3t_a(ji,jj,jk) = te3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 401 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 402 & ) * r1_e1e2t(ji,jj) 403 END DO 404 END DO 405 END DO 406 ! ! d - thickness diffusion transport: boundary conditions 407 ! (stored for tracer advction and continuity equation) 408 CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp) 409 410 ! 4 - Time stepping of baroclinic scale factors 411 ! --------------------------------------------- 412 ! Leapfrog time stepping 413 ! ~~~~~~~~~~~~~~~~~~~~~~ 414 CALL lbc_lnk( te3t_a(:,:,:), 'T', 1._wp ) 415 te3t_a(:,:,:) = te3t_b(:,:,:) + z2dt * tmask(:,:,:) * te3t_a(:,:,:) 416 417 ! Maximum deformation control 418 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 419 ze3t(:,:,jpk) = 0._wp 420 DO jk = 1, jpkm1 421 ze3t(:,:,jk) = te3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 422 END DO 423 z_tmax = MAXVAL( ze3t(:,:,:) ) 424 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 425 z_tmin = MINVAL( ze3t(:,:,:) ) 426 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 427 ! - ML - test: for the moment, stop simulation for too large e3_t variations 428 IF( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN 429 IF( lk_mpp ) THEN 430 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 431 CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 432 ELSE 433 ijk_max = MAXLOC( ze3t(:,:,:) ) 434 ijk_max(1) = ijk_max(1) + nimpp - 1 435 ijk_max(2) = ijk_max(2) + njmpp - 1 436 ijk_min = MINLOC( ze3t(:,:,:) ) 437 ijk_min(1) = ijk_min(1) + nimpp - 1 438 ijk_min(2) = ijk_min(2) + njmpp - 1 439 ENDIF 440 IF (lwp) THEN 441 WRITE(numout, *) 'MAX( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 442 WRITE(numout, *) 'at i, j, k=', ijk_max 443 WRITE(numout, *) 'MIN( te3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 444 WRITE(numout, *) 'at i, j, k=', ijk_min 445 CALL ctl_warn('MAX( ABS( te3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 446 ENDIF 447 ENDIF 448 ! - ML - end test 449 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 450 te3t_a(:,:,:) = MIN( te3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 451 te3t_a(:,:,:) = MAX( te3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 452 453 ! 454 ! "tilda" change in the after scale factor 455 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 456 DO jk = 1, jpkm1 457 dte3t_a(:,:,jk) = te3t_a(:,:,jk) - te3t_b(:,:,jk) 458 END DO 459 ! III - Barotropic repartition of the sea surface height over the baroclinic profile 460 ! ================================================================================== 461 ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 462 ! - ML - baroclinicity error should be better treated in the future 463 ! i.e. locally and not spread over the water column. 464 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 465 zht(:,:) = 0._wp 466 DO jk = 1, jpkm1 467 zht(:,:) = zht(:,:) + te3t_a(:,:,jk) * tmask(:,:,jk) 468 END DO 469 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 470 DO jk = 1, jpkm1 471 dte3t_a(:,:,jk) = dte3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 472 END DO 473 474 ENDIF 475 476 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 477 ! ! ---baroclinic part--------- ! 478 DO jk = 1, jpkm1 479 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dte3t_a(:,:,jk) * tmask(:,:,jk) 480 END DO 481 ENDIF 482 483 IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging 484 ! 485 IF( lwp ) WRITE(numout, *) 'kt =', kt 486 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 487 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 488 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 489 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(te3t_a))) =', z_tmax 490 END IF 491 ! 492 zht(:,:) = 0.0_wp 493 DO jk = 1, jpkm1 494 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 495 END DO 496 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 497 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 498 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 499 ! 500 zht(:,:) = 0.0_wp 501 DO jk = 1, jpkm1 502 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 503 END DO 504 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 505 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 506 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 507 ! 508 zht(:,:) = 0.0_wp 509 DO jk = 1, jpkm1 510 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 511 END DO 512 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 513 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 514 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 515 ! 516 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) 517 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 518 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 519 ! 520 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) ) 521 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 522 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 523 ! 524 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) ) 525 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 526 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 527 END IF 528 529 ! *********************************** ! 530 ! After scale factors at u- v- points ! 531 ! *********************************** ! 532 533 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 534 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 535 536 ! *********************************** ! 537 ! After depths at u- v points ! 538 ! *********************************** ! 539 540 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 541 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 542 DO jk = 2, jpkm1 543 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 544 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 545 END DO 546 ! ! Inverse of the local depth 547 !!gm BUG ? don't understand the use of umask_i here ..... 226 ! ! ------------------! 227 ! ! z_star coordinate ! (and barotropic z-tilde part) 228 ! ! ------------------! 229 ! 230 ! !== after ssh ==! (u- and v-points) 231 DO jj = 2, jpjm1 ; DO ji = 2, jpim1 232 zsshu_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji+1,jj) ) * ssumask(ji,jj) 233 zsshv_h(ji,jj) = 0.5_wp * ( ssha(ji,jj) + ssha(ji,jj+1) ) * ssvmask(ji,jj) 234 END DO ; END DO 235 CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp ) 236 ! 237 ! !== after depths and its inverse ==! 238 hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:) 239 hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:) 548 240 r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 549 241 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 242 ! 243 ! !== after scale factors ==! (e3t , e3u , e3v) 244 zssht_h(:,:) = ssha (:,:) * r1_ht_0(:,:) ! t-point 245 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! u-point 246 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! v-point 247 DO jk = 1, jpkm1 248 e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 249 e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 250 e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 251 END DO 550 252 ! 551 253 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 581 283 ! 582 284 INTEGER :: ji, jj, jk ! dummy loop indices 583 REAL(wp) :: zcoef, ze3f ! local scalar285 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h, zsshf_h 584 286 !!---------------------------------------------------------------------- 585 287 ! … … 594 296 ENDIF 595 297 ! 596 ! Time filter and swap of scale factors 597 ! ===================================== 598 ! - ML - e3(t/u/v)_b are allready computed in dynnxt. 599 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 600 IF( l_1st_euler ) THEN 601 te3t_n(:,:,:) = te3t_a(:,:,:) 602 ELSE 603 DO jk = 1, jpk 604 DO jj = 1, jpj 605 DO ji = 1, jpi 606 ze3f = te3t_n(ji,jj,jk) & 607 & + rn_atfp * ( te3t_b(ji,jj,jk) - 2.0_wp * te3t_n(ji,jj,jk) + te3t_a(ji,jj,jk) ) 608 te3t_b(ji,jj,jk) = ze3f 609 te3t_n(ji,jj,jk) = te3t_a(ji,jj,jk) 610 END DO 611 END DO 612 END DO 613 ENDIF 614 ENDIF 615 gdept_b(:,:,:) = gdept_n(:,:,:) 616 gdepw_b(:,:,:) = gdepw_n(:,:,:) 617 618 e3t_n(:,:,:) = e3t_a(:,:,:) 619 e3u_n(:,:,:) = e3u_a(:,:,:) 620 e3v_n(:,:,:) = e3v_a(:,:,:) 621 622 ! Compute all missing vertical scale factor and depths 623 ! ==================================================== 624 ! Horizontal scale factor interpolations 625 ! -------------------------------------- 298 ! Swap and Compute all missing vertical scale factor and depths 299 ! ============================================================= 626 300 ! - ML - e3u_b and e3v_b are allready computed in dynnxt 627 301 ! - JC - hu_b, hv_b, hur_b, hvr_b also 628 629 CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F' ) 630 631 ! Vertical scale factor interpolations 632 CALL dom_vvl_interpol( e3t_n(:,:,:), e3w_n(:,:,:), 'W' ) 633 CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' ) 634 CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' ) 635 CALL dom_vvl_interpol( e3t_b(:,:,:), e3w_b(:,:,:), 'W' ) 636 CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' ) 637 CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' ) 638 639 ! t- and w- points depth (set the isf depth as it is in the initial step) 640 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 641 gdepw_n(:,:,1) = 0.0_wp 642 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 643 DO jk = 2, jpk 644 DO jj = 1,jpj 645 DO ji = 1,jpi 646 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 647 ! 1 for jk = mikt 648 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 649 gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1) 650 gdept_n(ji,jj,jk) = zcoef * ( gdepw_n(ji,jj,jk ) + 0.5 * e3w_n(ji,jj,jk) ) & 651 & + (1-zcoef) * ( gdept_n(ji,jj,jk-1) + e3w_n(ji,jj,jk) ) 652 gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj) 653 END DO 654 END DO 655 END DO 656 302 ! 303 ! - GM - to be updated : e3f_n, e3w_n , e3uw_n , e3vw_n 304 ! e3w_b , e3uw_b , e3vw_b 305 ! gdept_n , gdepw_n , gde3w_n 306 ! ht_n 307 ! to be swap : hu_a , hv_a , r1_hu_a , r1_hv_a 308 ! 657 309 ! Local depth and Inverse of the local depth of the water 658 310 ! ------------------------------------------------------- 311 ! ! swap of depth and scale factors 312 ! ! =============================== 313 DO jk = 1, jpkm1 ! swap n--> a 314 gdept_b(:,:,jk) = gdept_n(:,:,jk) ! depth at t and w 315 gdepw_b(:,:,jk) = gdepw_n(:,:,jk) 316 e3t_n (:,:,jk) = e3t_a (:,:,jk) ! e3t, e3u, e3v 317 e3u_n (:,:,jk) = e3u_a (:,:,jk) 318 e3v_n (:,:,jk) = e3v_a (:,:,jk) 319 END DO 320 ht_n(:,:) = ht_0(:,:) + sshn(:,:) ! ocean thickness 321 ! 659 322 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 660 323 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 661 324 ! 662 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 663 DO jk = 2, jpkm1 664 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 325 ! !== before : 326 ! !* ssh at u- and v-points) 327 DO jj = 2, jpjm1 ; DO ji = 2, jpim1 328 zsshu_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) 329 zsshv_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) 330 END DO ; END DO 331 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) 332 ! 333 ! !* e3w_b , e3uw_b , e3vw_b 334 zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:) ! w-point 335 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! uw-point 336 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! vw-point 337 DO jk = 1, jpkm1 338 e3w_b(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 339 e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 340 e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 665 341 END DO 666 342 ! 343 ! !== now : 344 ! !* ssh at u- and v-points) 345 DO jj = 1, jpjm1 ; DO ji = 1, jpim1 ! start from 1 for f-point 346 zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) 347 zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) 348 zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 349 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 350 END DO ; END DO 351 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp ) 352 ! 353 ! !* e3w_n , e3uw_n , e3vw_n, e3f_n 354 zssht_h(:,:) = sshn (:,:) * r1_ht_0(:,:) ! t- & w-point 355 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) ! uw-point 356 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) ! vw-point 357 zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:) ! f-point 358 DO jk = 1, jpkm1 359 e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 360 e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) 361 e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) 362 e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) ) 363 END DO 364 ! 365 zssht_h(:,:) = 1._wp + sshn (:,:) * r1_ht_0(:,:) ! t-point 366 ! 367 IF( ln_isfcav ) THEN ! ISF cavities : ssh scaling not applied over the iceshelf thickness 368 DO jk = 1, jpkm1 369 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 370 gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 371 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn (:,:) 372 END DO 373 ELSE ! no ISF cavities 374 DO jk = 1, jpkm1 375 gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) 376 gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) 377 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 378 END DO 379 ENDIF 380 ! 667 381 ! write restart file 668 382 ! ================== … … 672 386 ! 673 387 END SUBROUTINE dom_vvl_sf_swp 674 675 676 SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )677 !!---------------------------------------------------------------------678 !! *** ROUTINE dom_vvl__interpol ***679 !!680 !! ** Purpose : interpolate scale factors from one grid point to another681 !!682 !! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0)683 !! - horizontal interpolation: grid cell surface averaging684 !! - vertical interpolation: simple averaging685 !!----------------------------------------------------------------------686 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pe3_in ! input e3 to be interpolated687 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe3_out ! output interpolated e3688 CHARACTER(LEN=*) , INTENT(in ) :: pout ! grid point of out scale factors689 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW'690 !691 INTEGER :: ji, jj, jk ! dummy loop indices692 REAL(wp) :: zlnwd ! =1./0. when ln_wd_il = T/F693 !!----------------------------------------------------------------------694 !695 IF(ln_wd_il) THEN696 zlnwd = 1.0_wp697 ELSE698 zlnwd = 0.0_wp699 END IF700 !701 SELECT CASE ( pout ) !== type of interpolation ==!702 !703 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean704 DO jk = 1, jpk705 DO jj = 1, jpjm1706 DO ji = 1, fs_jpim1 ! vector opt.707 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) &708 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) &709 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )710 END DO711 END DO712 END DO713 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )714 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)715 !716 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean717 DO jk = 1, jpk718 DO jj = 1, jpjm1719 DO ji = 1, fs_jpim1 ! vector opt.720 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) &721 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) &722 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )723 END DO724 END DO725 END DO726 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )727 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)728 !729 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean730 DO jk = 1, jpk731 DO jj = 1, jpjm1732 DO ji = 1, fs_jpim1 ! vector opt.733 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &734 & * r1_e1e2f(ji,jj) &735 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) &736 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )737 END DO738 END DO739 END DO740 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )741 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)742 !743 CASE( 'W' ) !* from T- to W-point : vertical simple mean744 !745 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)746 ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing747 !!gm BUG? use here wmask in case of ISF ? to be checked748 DO jk = 2, jpk749 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &750 & * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) &751 & + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &752 & * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) )753 END DO754 !755 CASE( 'UW' ) !* from U- to UW-point : vertical simple mean756 !757 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)758 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing759 !!gm BUG? use here wumask in case of ISF ? to be checked760 DO jk = 2, jpk761 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &762 & * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) &763 & + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &764 & * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) )765 END DO766 !767 CASE( 'VW' ) !* from V- to VW-point : vertical simple mean768 !769 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)770 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing771 !!gm BUG? use here wvmask in case of ISF ? to be checked772 DO jk = 2, jpk773 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) &774 & * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) &775 & + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) &776 & * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) )777 END DO778 END SELECT779 !780 END SUBROUTINE dom_vvl_interpol781 388 782 389 … … 804 411 IF( ln_rstart ) THEN !* Read the restart file 805 412 CALL rst_read_open ! open the restart file if necessary 806 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios )807 413 ! 808 id1 = iom_varid( numror, 'e3t_b' , ldstop = .FALSE. ) 809 id2 = iom_varid( numror, 'e3t_n' , ldstop = .FALSE. ) 810 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 811 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 812 id5 = iom_varid( numror, 'hdiv_lf' , ldstop = .FALSE. ) 813 ! ! --------- ! 814 ! ! all cases ! 815 ! ! --------- ! 414 id1 = iom_varid( numror, 'sshb' , ldstop = .FALSE. ) 415 id2 = iom_varid( numror, 'sshn' , ldstop = .FALSE. ) 416 ! 816 417 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 817 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 818 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 819 ! needed to restart if land processor not computed 820 IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files' 821 WHERE ( tmask(:,:,:) == 0.0_wp ) 822 e3t_n(:,:,:) = e3t_0(:,:,:) 823 e3t_b(:,:,:) = e3t_0(:,:,:) 824 END WHERE 418 IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files' 419 ! 420 !!gm Question: use jpdom_data above to read data over jpi x jpj (like is dom_hgr_read and dom_zgr_read) 421 !! so that it will work with land processor suppression 422 ! CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn, ldxios = lrxios ) 423 ! CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb, ldxios = lrxios ) 424 !!gm 425 CALL iom_get( numror, jpdom_data, 'sshn' , sshn, ldxios = lrxios ) 426 CALL iom_get( numror, jpdom_data, 'sshb' , sshb, ldxios = lrxios ) 427 !!gm end 825 428 IF( l_1st_euler ) THEN 826 e3t_b(:,:,:) = e3t_n(:,:,:)429 sshb(:,:) = sshn(:,:) 827 430 ENDIF 828 431 ELSE IF( id1 > 0 ) THEN 829 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'830 IF(lwp) write(numout,*) ' e3t_n set equal to e3t_b.'831 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 832 CALL iom_get( numror, jpdom_ autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )833 e3t_n(:,:,:) = e3t_b(:,:,:)432 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files' 433 IF(lwp) write(numout,*) ' set sshn = sshb and force l_1st_euler = true' 434 !!gm CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios ) 435 CALL iom_get( numror, jpdom_data, 'sshb', sshb, ldxios = lrxios ) 436 sshn(:,:) = sshb(:,:) 834 437 l_1st_euler = .TRUE. 835 438 ELSE IF( id2 > 0 ) THEN 836 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files' 837 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 838 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 839 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios ) 840 e3t_b(:,:,:) = e3t_n(:,:,:) 439 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files' 440 IF(lwp) write(numout,*) 'set sshb = sshn and force l_1st_euler = true' 441 CALL iom_get( numror, jpdom_data, 'sshn', sshb, ldxios = lrxios ) 442 sshb(:,:) = sshn(:,:) 841 443 l_1st_euler = .TRUE. 842 444 ELSE 843 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file' 844 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 845 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 846 DO jk = 1, jpk 847 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 848 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 849 & + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk)) 850 END DO 851 e3t_b(:,:,:) = e3t_n(:,:,:) 445 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file' 446 IF(lwp) write(numout,*) 'set sshb = sshn = 0 and force l_1st_euler = true' 447 sshb(:,:) = 0._wp 448 sshn(:,:) = 0._wp 852 449 l_1st_euler = .TRUE. 853 450 ENDIF 854 ! ! ----------- !855 IF( ln_vvl_zstar ) THEN ! z_star case !856 ! ! ----------- !857 IF( MIN( id3, id4 ) > 0 ) THEN858 CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )859 ENDIF860 ! ! ----------------------- !861 ELSE ! z_tilde and layer cases !862 ! ! ----------------------- !863 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist864 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lrxios )865 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lrxios )866 ELSE ! one at least array is missing867 te3t_b(:,:,:) = 0.0_wp868 te3t_n(:,:,:) = 0.0_wp869 ENDIF870 ! ! ------------ !871 IF( ln_vvl_ztilde ) THEN ! z_tilde case !872 ! ! ------------ !873 IF( id5 > 0 ) THEN ! required array exists874 CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )875 ELSE ! array is missing876 hdiv_lf(:,:,:) = 0.0_wp877 ENDIF878 ENDIF879 ENDIF880 !881 451 ELSE !* Initialize at "rest" 882 452 ! 883 884 453 IF( ll_wd ) THEN ! MJB ll_wd edits start here - these are essential 885 454 ! 886 IF( cn_cfg == 'wad' ) THEN 887 ! Wetting and drying test case 455 IF( cn_cfg == 'wad' ) THEN ! Wetting and drying test case 888 456 CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb ) 889 tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones457 tsn (:,:,:,:) = tsb (:,:,:,:) ! set now values from to before ones 890 458 sshn (:,:) = sshb(:,:) 891 459 un (:,:,:) = ub (:,:,:) 892 460 vn (:,:,:) = vb (:,:,:) 893 ELSE 894 ! if not test case 461 ELSE ! Not the test case 895 462 sshn(:,:) = -ssh_ref 896 463 sshb(:,:) = -ssh_ref 897 464 ! 898 465 DO jj = 1, jpj 899 466 DO ji = 1, jpi 900 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 901 467 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 902 468 sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 903 469 sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 904 470 ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) ) 905 471 ENDIF 906 ENDDO 907 ENDDO 908 ENDIF !If test case else 909 910 ! Adjust vertical metrics for all wad 911 DO jk = 1, jpk 912 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) & 913 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 914 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) 915 END DO 916 e3t_b(:,:,:) = e3t_n(:,:,:) 917 472 END DO 473 END DO 474 ENDIF ! If test case else 475 ! 918 476 DO ji = 1, jpi 919 477 DO jj = 1, jpj 920 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN478 IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN 921 479 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 922 480 ENDIF … … 926 484 ELSE 927 485 ! 928 ! Just to read set ssh in fact, called latter once vertical grid 929 ! is set up: 486 ! Just to read set ssh in fact, called latter once vertical grid is set up: 930 487 ! CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb ) 931 ! ! 932 ! DO jk=1,jpk 933 ! e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) & 934 ! & / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk) 935 ! END DO 936 ! e3t_n(:,:,:) = e3t_b(:,:,:) 937 sshn(:,:)=0._wp 938 e3t_n(:,:,:)=e3t_0(:,:,:) 939 e3t_b(:,:,:)=e3t_0(:,:,:) 488 sshn(:,:) = 0._wp 489 sshb(:,:) = 0._wp 940 490 ! 941 END IF ! end of ll_wd edits942 943 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN944 te3t_b(:,:,:) = 0._wp945 te3t_n(:,:,:) = 0._wp946 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp947 491 END IF 492 ! 948 493 ENDIF 949 494 ! 950 495 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 951 496 ! ! =================== 952 IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 953 IF( lwxios ) CALL iom_swap( cwxios_context ) 954 ! ! --------- ! 955 ! ! all cases ! 956 ! ! --------- ! 957 CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios ) 958 CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios ) 959 ! ! ----------------------- ! 960 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! 961 ! ! ----------------------- ! 962 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', te3t_b(:,:,:), ldxios = lwxios) 963 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', te3t_n(:,:,:), ldxios = lwxios) 964 END IF 965 ! ! -------------! 966 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 967 ! ! ------------ ! 968 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios) 969 ENDIF 970 ! 971 IF( lwxios ) CALL iom_swap( cxios_context ) 497 498 !!gm DO NOTHING, sshb and sshn are written in restart.F90 499 972 500 ENDIF 973 501 ! … … 1024 552 ENDIF 1025 553 ! 554 555 !!gm 556 IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar ) CALL ctl_stop( 'z-tilde NOT available in this branch' ) 557 !!gm 558 559 ! 1026 560 ioptio = 0 ! Parameter control 1027 561 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true. … … 1047 581 END SUBROUTINE dom_vvl_ctl 1048 582 583 584 SUBROUTINE ssh2e3_now 585 !!---------------------------------------------------------------------- 586 !! *** ROUTINE ssh2e3_now *** 587 !!---------------------------------------------------------------------- 588 INTEGER :: ji, jj, jk 589 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h, zsshf_h 590 !!---------------------------------------------------------------------- 591 ! 592 ! !== ssh at u- and v-points ==! 593 ! 594 DO jj = 1, jpjm1 ! start from 1 due to f-point 595 DO ji = 1, jpim1 596 zsshu_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji+1,jj ) ) * ssumask(ji,jj) 597 zsshv_h(ji,jj) = 0.50_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) ) * ssvmask(ji,jj) 598 zsshf_h(ji,jj) = 0.25_wp * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 599 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) * ssfmask(ji,jj) 600 END DO 601 END DO 602 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp ) 603 ! 604 ! !== ht, hu and hv == ! (and their inverse) 605 ! 606 ht_n (:,:) = ht_0(:,:) + sshn (:,:) 607 hu_n (:,:) = hu_0(:,:) + zsshu_h(:,:) 608 hv_n (:,:) = hv_0(:,:) + zsshv_h(:,:) 609 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) ! ss mask mask due to ISF 610 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 611 ! 612 ! !== ssh / h factor at t-, u- ,v- & f-points ==! 613 ! 614 zssht_h(:,:) = sshn (:,:) * r1_ht_0(:,:) 615 zsshu_h(:,:) = zsshu_h(:,:) * r1_hu_0(:,:) 616 zsshv_h(:,:) = zsshv_h(:,:) * r1_hv_0(:,:) 617 zsshf_h(:,:) = zsshf_h(:,:) * r1_hf_0(:,:) 618 ! 619 ! !== e3t, e3w , e3u, e3uw , e3v, e3vw , and e3f ==! 620 ! 621 DO jk = 1, jpkm1 622 e3t_n(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 623 e3w_n(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 624 ! 625 e3u_n(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) ) 626 e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) ) 627 ! 628 e3v_n(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) ) 629 e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) ) 630 ! 631 e3f_n(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) ) 632 END DO 633 ! 634 ! !== depth of t- and w-points ==! 635 ! 636 zssht_h(:,:) = 1._wp + zssht_h(:,:) ! = 1 + sshn / ht_0 637 ! 638 IF( ln_isfcav ) THEN ! ISF cavities : ssh scaling not applied over the iceshelf thickness 639 DO jk = 1, jpkm1 640 gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 641 gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:) 642 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 643 END DO 644 ELSE ! no ISF cavities 645 DO jk = 1, jpkm1 646 gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:) 647 gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:) 648 gde3w_n(:,:,jk) = gdept_n(:,:,jk) - sshn(:,:) 649 END DO 650 ENDIF 651 ! 652 END SUBROUTINE ssh2e3_now 653 654 655 SUBROUTINE ssh2e3_before 656 !!---------------------------------------------------------------------- 657 !! *** ROUTINE ssh2e3_before *** 658 !!---------------------------------------------------------------------- 659 INTEGER :: ji, jj, jk 660 REAL(wp), DIMENSION(jpi,jpj) :: zssht_h, zsshu_h, zsshv_h 661 !!---------------------------------------------------------------------- 662 ! 663 ! !== ssh at u- and v-points ==! 664 DO jj = 2, jpjm1 665 DO ji = 2, jpim1 666 zsshu_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji+1,jj ) ) * ssumask(ji,jj) 667 zsshv_h(ji,jj) = 0.5_wp * ( sshb(ji ,jj) + sshb(ji ,jj+1) ) * ssvmask(ji,jj) 668 END DO 669 END DO 670 CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp ) 671 ! 672 ! !== ht, hu and hv == ! (and their inverse) 673 hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:) 674 hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:) 675 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) 676 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 677 ! 678 ! 679 ! !== ssh / h factor at t-, u- ,v- & f-points ==! 680 zssht_h(:,:) = sshb (:,:) * r1_ht_0(:,:) 681 zsshu_h (:,:) = zsshu_h(:,:) * r1_hu_0(:,:) 682 zsshv_h (:,:) = zsshv_h(:,:) * r1_hv_0(:,:) 683 ! 684 ! !== e3t, e3w , e3u, e3uw , and e3v, e3vw ==! 685 DO jk = 1, jpkm1 686 e3t_b(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 687 e3w_b(:,:,jk) = e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) ) 688 ! 689 e3u_b(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h (:,:) * umask(:,:,jk) ) 690 e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h (:,:) * wumask(:,:,jk) ) 691 ! 692 e3v_b(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h (:,:) * vmask(:,:,jk) ) 693 e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h (:,:) * wvmask(:,:,jk) ) 694 END DO 695 ! 696 END SUBROUTINE ssh2e3_before 697 1049 698 !!====================================================================== 1050 699 END MODULE domvvl
Note: See TracChangeset
for help on using the changeset viewer.
