Changeset 12740 for NEMO/trunk/tests/VORTEX/MY_SRC/domvvl.F90
- Timestamp:
- 2020-04-12T11:03:06+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/tests/VORTEX/MY_SRC/domvvl.F90
r12489 r12740 63 63 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 64 64 65 !! * Substitutions 66 # include "do_loop_substitute.h90" 65 67 !!---------------------------------------------------------------------- 66 68 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 188 190 gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb) 189 191 gdepw(:,:,1,Kbb) = 0.0_wp 190 DO jk = 2, jpk ! vertical sum 191 DO jj = 1,jpj 192 DO ji = 1,jpi 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 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 196 196 !!gm ??????? BUG ? gdept(:,:,:,Kmm) as well as gde3w does not include the thickness of ISF ?? 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 DO 206 END DO 207 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 208 206 ! 209 207 ! !== thickness of the water column !! (ocean portion only) … … 240 238 ENDIF 241 239 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 242 DO jj = 1, jpj 243 DO ji = 1, jpi 240 DO_2D_11_11 244 241 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 245 IF( ABS(gphit(ji,jj)) >= 6.) THEN 246 ! values outside the equatorial band and transition zone (ztilde) 247 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 248 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 249 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 250 ! values inside the equatorial band (ztilde as zstar) 251 frq_rst_e3t(ji,jj) = 0.0_wp 252 frq_rst_hdv(ji,jj) = 1.0_wp / rn_Dt 253 ELSE ! transition band (2.5 to 6 degrees N/S) 254 ! ! (linearly transition from z-tilde to z-star) 255 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 256 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 257 & * 180._wp / 3.5_wp ) ) 258 frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt) & 259 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp & 260 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 261 & * 180._wp / 3.5_wp ) ) 262 ENDIF 263 END DO 264 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 265 261 IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN 266 262 IF( nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2 … … 357 353 END DO 358 354 ! 359 IF( ln_vvl_ztilde .OR. ln_vvl_layer.AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate !360 ! ! ------baroclinic part------ !355 IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 356 ! ! ------baroclinic part------ ! 361 357 ! I - initialization 362 358 ! ================== … … 411 407 zwu(:,:) = 0._wp 412 408 zwv(:,:) = 0._wp 413 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 414 DO jj = 1, jpjm1 415 DO ji = 1, jpim1 ! vector opt. 416 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 417 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 418 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 419 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 420 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 421 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 422 END DO 423 END DO 424 END DO 425 DO jj = 1, jpj ! b - correction for last oceanic u-v points 426 DO ji = 1, jpi 427 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 428 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 429 END DO 430 END DO 431 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 432 DO jj = 2, jpjm1 433 DO ji = 2, jpim1 ! vector opt. 434 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 435 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 436 & ) * r1_e1e2t(ji,jj) 437 END DO 438 END DO 439 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 440 426 ! ! d - thickness diffusion transport: boundary conditions 441 427 ! (stored for tracer advction and continuity equation) … … 444 430 ! 4 - Time stepping of baroclinic scale factors 445 431 ! --------------------------------------------- 446 ! Leapfrog time stepping447 ! ~~~~~~~~~~~~~~~~~~~~~~448 432 CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp ) 449 433 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:) … … 646 630 ! Horizontal scale factor interpolations 647 631 ! -------------------------------------- 648 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are al lready computed in dynnxt632 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 649 633 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 650 634 … … 663 647 gdepw(:,:,1,Kmm) = 0.0_wp 664 648 gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm) 665 DO jk = 2, jpk 666 DO jj = 1,jpj 667 DO ji = 1,jpi 668 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 669 ! 1 for jk = mikt 670 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 671 gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm) 672 gdept(ji,jj,jk,Kmm) = zcoef * ( gdepw(ji,jj,jk ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) ) & 673 & + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) + e3w(ji,jj,jk,Kmm) ) 674 gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm) 675 END DO 676 END DO 677 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 678 658 679 659 ! Local depth and Inverse of the local depth of the water … … 722 702 ! 723 703 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 724 DO jk = 1, jpk 725 DO jj = 1, jpjm1 726 DO ji = 1, jpim1 ! vector opt. 727 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj) & 728 & * ( e1e2t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 729 & + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 730 END DO 731 END DO 732 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 733 709 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp ) 734 710 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 735 711 ! 736 712 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 737 DO jk = 1, jpk 738 DO jj = 1, jpjm1 739 DO ji = 1, jpim1 ! vector opt. 740 pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj) & 741 & * ( e1e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 742 & + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 743 END DO 744 END DO 745 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 746 718 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp ) 747 719 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 748 720 ! 749 721 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 750 DO jk = 1, jpk 751 DO jj = 1, jpjm1 752 DO ji = 1, jpim1 ! vector opt. 753 pe3_out(ji,jj,jk) = 0.5_wp * ( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) & 754 & * r1_e1e2f(ji,jj) & 755 & * ( e1e2u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 756 & + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 757 END DO 758 END DO 759 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 760 728 CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp ) 761 729 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) … … 832 800 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 833 801 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 802 ! 834 803 ! ! --------- ! 835 804 ! ! all cases ! 836 805 ! ! --------- ! 806 ! 837 807 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 838 808 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) … … 850 820 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files' 851 821 IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.' 852 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'822 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 853 823 CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios ) 854 824 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) … … 857 827 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files' 858 828 IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.' 859 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'829 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 860 830 CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios ) 861 831 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) … … 864 834 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file' 865 835 IF(lwp) write(numout,*) 'Compute scale factor from sshn' 866 IF(lwp) write(numout,*) 'l_1st_euler is forced to .true.'836 IF(lwp) write(numout,*) 'l_1st_euler is forced to true' 867 837 DO jk = 1, jpk 868 838 e3t(:,:,jk,Kmm) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) & … … 917 887 ssh(:,:,Kbb) = -ssh_ref 918 888 919 DO jj = 1, jpj 920 DO ji = 1, jpi 921 IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth 922 ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) ) 923 ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) ) 924 ENDIF 925 ENDDO 926 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 927 895 ENDIF !If test case else 928 896 … … 935 903 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm) 936 904 937 DO ji = 1, jpi 938 DO jj = 1, jpj 939 IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN 940 CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' ) 941 ENDIF 942 END DO 943 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 944 910 ! 945 911 ELSE
Note: See TracChangeset
for help on using the changeset viewer.