Changeset 13197 for NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/domvvl.F90
- Timestamp:
- 2020-07-01T16:09:00+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/tests/CANAL/MY_SRC/domvvl.F90
r12489 r13197 37 37 38 38 PUBLIC dom_vvl_init ! called by domain.F90 39 PUBLIC dom_vvl_zgr ! called by isfcpl.F90 39 40 PUBLIC dom_vvl_sf_nxt ! called by step.F90 40 41 PUBLIC dom_vvl_sf_update ! called by step.F90 … … 62 63 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 63 64 65 !! * Substitutions 66 # include "do_loop_substitute.h90" 64 67 !!---------------------------------------------------------------------- 65 68 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 116 119 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 117 120 ! 121 IF(lwp) WRITE(numout,*) 122 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 123 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 124 ! 125 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer) 126 ! 127 ! ! Allocate module arrays 128 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 129 ! 130 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 131 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' ) 132 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 133 ! 134 CALL dom_vvl_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column 135 ! 136 END SUBROUTINE dom_vvl_init 137 ! 138 SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa) 139 !!---------------------------------------------------------------------- 140 !! *** ROUTINE dom_vvl_init *** 141 !! 142 !! ** Purpose : Interpolation of all scale factors, 143 !! depths and water column heights 144 !! 145 !! ** Method : - interpolate scale factors 146 !! 147 !! ** Action : - e3t_(n/b) and tilde_e3t_(n/b) 148 !! - Regrid: e3(u/v)_n 149 !! e3(u/v)_b 150 !! e3w_n 151 !! e3(u/v)w_b 152 !! e3(u/v)w_n 153 !! gdept_n, gdepw_n and gde3w_n 154 !! - h(t/u/v)_0 155 !! - frq_rst_e3t and frq_rst_hdv 156 !! 157 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 158 !!---------------------------------------------------------------------- 159 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa 160 !!---------------------------------------------------------------------- 118 161 INTEGER :: ji, jj, jk 119 162 INTEGER :: ii0, ii1, ij0, ij1 120 163 REAL(wp):: zcoef 121 164 !!---------------------------------------------------------------------- 122 !123 IF(lwp) WRITE(numout,*)124 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'125 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'126 !127 CALL dom_vvl_ctl ! choose vertical coordinate (z_star, z_tilde or layer)128 !129 ! ! Allocate module arrays130 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )131 !132 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf133 CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )134 e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all135 165 ! 136 166 ! !== Set of all other vertical scale factors ==! (now and before) … … 160 190 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 161 191 gdepw(:,:,1,Kbb) = 0.0_wp 162 DO jk = 2, jpk ! vertical sum 163 DO jj = 1,jpj 164 DO ji = 1,jpi 165 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 166 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 167 ! ! 0.5 where jk = mikt 192 DO_3D_11_11( 2, jpk ) 193 ! zcoef = tmask - wmask ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 194 ! ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf) 195 ! ! 0.5 where jk = mikt 168 196 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 169 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 170 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 171 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 172 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 173 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 174 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 175 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 176 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 177 END DO 178 END DO 179 END DO 197 zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) ) 198 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 199 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm)) & 200 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm)) 201 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 202 gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb) 203 gdept(ji,jj,jk,Kbb) = zcoef * ( gdepw(ji,jj,jk ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb)) & 204 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) + e3w(ji,jj,jk,Kbb)) 205 END_3D 180 206 ! 181 207 ! !== thickness of the water column !! (ocean portion only) … … 212 238 ENDIF 213 239 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 214 DO jj = 1, jpj 215 DO ji = 1, jpi 240 DO_2D_11_11 216 241 !!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.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 220 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_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.0_wp 224 frq_rst_hdv(ji,jj) = 1.0_wp / rn_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.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 228 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 229 & * 180._wp / 3.5_wp ) ) 230 frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) & 231 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp & 232 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 233 & * 180._wp / 3.5_wp ) ) 234 ENDIF 235 END DO 236 END DO 242 IF( ABS(gphit(ji,jj)) >= 6.) THEN 243 ! values outside the equatorial band and transition zone (ztilde) 244 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 245 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 246 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 247 ! values inside the equatorial band (ztilde as zstar) 248 frq_rst_e3t(ji,jj) = 0.0_wp 249 frq_rst_hdv(ji,jj) = 1.0_wp / rn_Dt 250 ELSE ! transition band (2.5 to 6 degrees N/S) 251 ! ! (linearly transition from z-tilde to z-star) 252 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 253 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 254 & * 180._wp / 3.5_wp ) ) 255 frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) & 256 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp & 257 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 258 & * 180._wp / 3.5_wp ) ) 259 ENDIF 260 END_2D 237 261 IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 238 262 IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 … … 264 288 ENDIF 265 289 ! 266 END SUBROUTINE dom_vvl_ init290 END SUBROUTINE dom_vvl_zgr 267 291 268 292 … … 329 353 END DO 330 354 ! 331 IF( ln_vvl_ztilde .OR. ln_vvl_layer.AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !332 ! ! ------baroclinic part------ !355 IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 356 ! ! ------baroclinic part------ ! 333 357 ! I - initialization 334 358 ! ================== … … 383 407 zwu(:,:) = 0._wp 384 408 zwv(:,:) = 0._wp 385 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 386 DO jj = 1, jpjm1 387 DO ji = 1, fs_jpim1 ! vector opt. 388 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 389 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 390 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 391 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 392 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 393 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 394 END DO 395 END DO 396 END DO 397 DO jj = 1, jpj ! b - correction for last oceanic u-v points 398 DO ji = 1, jpi 399 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 400 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 401 END DO 402 END DO 403 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 404 DO jj = 2, jpjm1 405 DO ji = fs_2, fs_jpim1 ! vector opt. 406 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 407 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 408 & ) * r1_e1e2t(ji,jj) 409 END DO 410 END DO 411 END DO 409 DO_3D_10_10( 1, jpkm1 ) 410 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 411 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 412 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 413 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 414 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 415 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 416 END_3D 417 DO_2D_11_11 418 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 419 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 420 END_2D 421 DO_3D_00_00( 1, jpkm1 ) 422 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 423 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 424 & ) * r1_e1e2t(ji,jj) 425 END_3D 412 426 ! ! d - thickness diffusion transport: boundary conditions 413 427 ! (stored for tracer advction and continuity equation) … … 416 430 ! 4 - Time stepping of baroclinic scale factors 417 431 ! --------------------------------------------- 418 ! Leapfrog time stepping419 ! ~~~~~~~~~~~~~~~~~~~~~~420 432 CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 421 433 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) … … 613 625 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 614 626 ENDIF 615 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)616 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)617 618 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)619 e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)620 e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)621 627 622 628 ! Compute all missing vertical scale factor and depths … … 641 647 gdepw(:,:,1,Kmm) = 0.0_wp 642 648 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 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(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 650 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 651 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 652 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 653 END DO 654 END DO 655 END DO 649 DO_3D_11_11( 2, jpk ) 650 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 651 ! 1 for jk = mikt 652 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 653 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 654 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 655 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 656 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 657 END_3D 656 658 657 659 ! Local depth and Inverse of the local depth of the water … … 700 702 ! 701 703 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 702 DO jk = 1, jpk 703 DO jj = 1, jpjm1 704 DO ji = 1, fs_jpim1 ! vector opt. 705 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 706 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 707 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 708 END DO 709 END DO 710 END DO 704 DO_3D_10_10( 1, jpk ) 705 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 706 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 707 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 708 END_3D 711 709 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 712 710 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 713 711 ! 714 712 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 715 DO jk = 1, jpk 716 DO jj = 1, jpjm1 717 DO ji = 1, fs_jpim1 ! vector opt. 718 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 719 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 720 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 721 END DO 722 END DO 723 END DO 713 DO_3D_10_10( 1, jpk ) 714 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 715 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 716 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 717 END_3D 724 718 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 725 719 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 726 720 ! 727 721 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 728 DO jk = 1, jpk 729 DO jj = 1, jpjm1 730 DO ji = 1, fs_jpim1 ! vector opt. 731 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 732 & * r1_e1e2f(ji,jj) & 733 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 734 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 735 END DO 736 END DO 737 END DO 722 DO_3D_10_10( 1, jpk ) 723 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 724 & * r1_e1e2f(ji,jj) & 725 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 726 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 727 END_3D 738 728 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 739 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) … … 810 800 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 811 801 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 802 ! 812 803 ! ! --------- ! 813 804 ! ! all cases ! 814 805 ! ! --------- ! 806 ! 815 807 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 816 808 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) … … 828 820 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 829 821 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 830 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'822 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 831 823 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 832 824 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) … … 835 827 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 836 828 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 837 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'829 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 838 830 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 839 831 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) … … 842 834 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 843 835 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 844 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'836 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 845 837 DO jk = 1, jpk 846 838 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & … … 895 887 ssh(:,:,Kbb) = -ssh_ref 896 888 897 DO jj = 1, jpj 898 DO ji = 1, jpi 899 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 900 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 901 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 902 ENDIF 903 ENDDO 904 ENDDO 889 DO_2D_11_11 890 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 891 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 892 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 893 ENDIF 894 END_2D 905 895 ENDIF !If test case else 906 896 … … 913 903 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 914 904 915 DO ji = 1, jpi 916 DO jj = 1, jpj 917 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 918 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 919 ENDIF 920 END DO 921 END DO 905 DO_2D_11_11 906 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 907 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 908 ENDIF 909 END_2D 922 910 ! 923 911 ELSE 924 912 ! 925 ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm) 926 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 927 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 913 ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 914 ! 915 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 916 ! 917 ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 928 918 ! 929 919 DO jk=1,jpk 930 e3t(:,:,jk,K mm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &931 & 932 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t_b!= 0 on land points920 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 921 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 922 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t(:,:,:,Kbb) != 0 on land points 933 923 END DO 934 924 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 935 ssh(:,: ,Kmm) = ssh(:,: ,Kbb)! needed later for gde3w925 ssh(:,:,Kmm) = ssh(:,:,Kbb) ! needed later for gde3w 936 926 ! 937 927 END IF ! end of ll_wd edits … … 1025 1015 ! 1026 1016 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 1027 IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )1028 1017 ! 1029 1018 IF(lwp) THEN ! Print the choice
Note: See TracChangeset
for help on using the changeset viewer.
