Changeset 10094
- Timestamp:
- 2018-09-06T15:24:19+02:00 (6 years ago)
- Location:
- NEMO/branches/UKMO/dev_r10037_dynvor_EEUV
- Files:
-
- 3 edited
-
cfgs/SHARED/namelist_ref (modified) (1 diff)
-
src/OCE/DYN/dynspg_ts.F90 (modified) (5 diffs)
-
src/OCE/DYN/dynvor.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/cfgs/SHARED/namelist_ref
r9950 r10094 868 868 ln_dynvor_enT = .false. ! energy conserving scheme (T-point) 869 869 ln_dynvor_eeT = .false. ! energy conserving scheme (een using e3t) 870 ln_dynvor_eeUV = .false.! energy conserving scheme (een using e3u and e3v) 870 871 ln_dynvor_een = .false. ! energy & enstrophy scheme 871 872 nn_een_e3f = 1 ! =0 e3f = mi(mj(e3t))/4 -
NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynspg_ts.F90
r9950 r10094 104 104 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 105 105 ! 106 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) &106 IF( ln_dynvor_een .OR. ln_dynvor_eeT .OR. ln_dynvor_eeUV) & 107 107 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 108 108 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) … … 156 156 REAL(wp) :: za0, za1, za2, za3 ! - - 157 157 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 158 REAL(wp) :: zmask_ne,zmask_nw,zmask_se,zmask_sw 159 REAL(wp) :: z1_huv_ne,z1_huv_nw,z1_huv_se,z1_huv_sw 158 160 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 159 161 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc … … 283 285 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 284 286 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 287 END DO 288 END DO 289 CASE( np_EEUV ) != EEN scheme using e3u and e3v (energy conserving scheme) 290 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 291 DO jj = 2, jpj 292 DO ji = 2, jpi 293 zmask_ne = MAX(ssumask(ji ,jj),ssvmask(ji,jj )) 294 zmask_nw = MAX(ssumask(ji-1,jj),ssvmask(ji,jj )) 295 zmask_se = MAX(ssumask(ji ,jj),ssvmask(ji,jj-1)) 296 zmask_sw = MAX(ssumask(ji-1,jj),ssvmask(ji,jj-1)) 297 z1_huv_ne = 2._wp * zmask_ne / ( hu_n(ji ,jj) + hv_n(ji,jj ) + 1._wp - zmask_ne ) 298 z1_huv_nw = 2._wp * zmask_nw / ( hu_n(ji-1,jj) + hv_n(ji,jj ) + 1._wp - zmask_nw ) 299 z1_huv_se = 2._wp * zmask_se / ( hu_n(ji ,jj) + hv_n(ji,jj-1) + 1._wp - zmask_se ) 300 z1_huv_sw = 2._wp * zmask_sw / ( hu_n(ji-1,jj) + hv_n(ji,jj-1) + 1._wp - zmask_sw ) 301 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_huv_ne 302 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_huv_nw 303 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_huv_se 304 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_huv_sw 285 305 END DO 286 306 END DO … … 426 446 END DO 427 447 ! 428 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)448 CASE( np_EET , np_EEN, np_EEUV ) ! energy & enstrophy scheme (using e3t or e3f or e3u and e3v) 429 449 DO jj = 2, jpjm1 430 450 DO ji = fs_2, fs_jpim1 ! vector opt. … … 1035 1055 END DO 1036 1056 ! 1037 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f)1057 CASE( np_EET , np_EEN, np_EEUV ) ! energy & enstrophy scheme (using e3t or e3f or e3u and e3v) 1038 1058 DO jj = 2, jpjm1 1039 1059 DO ji = fs_2, fs_jpim1 ! vector opt. -
NEMO/branches/UKMO/dev_r10037_dynvor_EEUV/src/OCE/DYN/dynvor.F90
r9950 r10094 56 56 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 57 57 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 58 LOGICAL, PUBLIC :: ln_dynvor_eeUV !: uv-point energy conserving scheme (EEUV) 58 59 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 59 60 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) … … 67 68 INTEGER, PUBLIC, PARAMETER :: np_ENT = 2 ! ENT scheme (t-point vorticity) 68 69 INTEGER, PUBLIC, PARAMETER :: np_EET = 3 ! EET scheme (EEN using e3t) 69 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 70 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 70 INTEGER, PUBLIC, PARAMETER :: np_EEUV = 4 ! EET scheme (EEN using e3u and e3v) 71 INTEGER, PUBLIC, PARAMETER :: np_EEN = 5 ! EEN scheme 72 INTEGER, PUBLIC, PARAMETER :: np_MIX = 6 ! MIX scheme 71 73 72 74 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity … … 84 86 85 87 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 88 REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6 86 89 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 87 90 REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12 … … 128 131 CASE( np_EET ) ; CALL vor_eeT( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 129 132 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 133 CASE( np_EEUV ) ; CALL vor_eeUV( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (een with e3u and e3v) 134 IF( ln_stcor ) CALL vor_eeUV( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 130 135 CASE( np_EEN ) ; CALL vor_een( kt, ncor, un , vn , ua, va ) ! energy & enstrophy scheme 131 136 IF( ln_stcor ) CALL vor_een( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend … … 141 146 CASE( np_ENT ) ; CALL vor_enT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (T-pts) 142 147 CASE( np_EET ) ; CALL vor_eeT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 148 CASE( np_EEUV ) ; CALL vor_eeUV( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (een with e3u and e3v) 143 149 CASE( np_ENE ) ; CALL vor_ene( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme 144 150 CASE( np_ENS, np_MIX ) ; CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! enstrophy conserving scheme … … 161 167 CALL vor_eeT( kt, ntot, un , vn , ua, va ) ! total vorticity trend 162 168 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 169 CASE( np_EEUV ) !* energy conserving scheme (een scheme using e3u and e3v) 170 CALL vor_eeUV( kt, ntot, un , vn , ua, va ) ! total vorticity trend 171 IF( ln_stcor ) CALL vor_eeUV( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 163 172 CASE( np_ENE ) !* energy conserving scheme 164 173 CALL vor_ene( kt, ntot, un , vn , ua, va ) ! total vorticity trend … … 806 815 807 816 817 SUBROUTINE vor_eeUV( kt, kvor, pun, pvn, pua, pva ) 818 !!---------------------------------------------------------------------- 819 !! *** ROUTINE vor_eeUV *** 820 !! 821 !! ** Purpose : Compute the now total vorticity trend and add it to 822 !! the general trend of the momentum equation. 823 !! 824 !! ** Method : Trend evaluated using now fields (centered in time) 825 !! and the Arakawa and Lamb (1980) flux form formulation : conserves 826 !! both the horizontal kinetic energy and the potential enstrophy 827 !! when horizontal divergence is zero (see the NEMO documentation) 828 !! Add this trend to the general momentum trend (ua,va). 829 !! 830 !! Modified version of EEN as suggested by Mike Bell. 831 !! 832 !! ** Action : - Update (ua,va) with the now vorticity term trend 833 !! 834 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 835 !! Bell, Johnson and Marshall, unpublished notes 2017 836 !!---------------------------------------------------------------------- 837 INTEGER , INTENT(in ) :: kt ! ocean time-step index 838 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 839 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun, pvn ! now velocities 840 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! total v-trend 841 ! 842 INTEGER :: ji, jj, jk ! dummy loop indices 843 INTEGER :: ierr ! local integer 844 REAL(wp) :: zua, zva ! local scalars 845 REAL(wp) :: zmsk, ze3f ! local scalars 846 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , zwz 847 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 848 !!---------------------------------------------------------------------- 849 ! 850 IF( kt == nit000 ) THEN 851 IF(lwp) WRITE(numout,*) 852 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 853 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 854 ENDIF 855 ! 856 ! ! =============== 857 DO jk = 1, jpkm1 ! Horizontal slab 858 ! ! =============== 859 ! 860 SELECT CASE( kvor ) !== vorticity considered ==! 861 CASE ( np_COR ) !* Coriolis (planetary vorticity) 862 DO jj = 1, jpjm1 863 DO ji = 1, fs_jpim1 ! vector opt. 864 zwz(ji,jj) = ff_f(ji,jj) 865 END DO 866 END DO 867 CASE ( np_RVO ) !* relative vorticity 868 DO jj = 1, jpjm1 869 DO ji = 1, fs_jpim1 ! vector opt. 870 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 871 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 872 END DO 873 END DO 874 CASE ( np_MET ) !* metric term 875 DO jj = 1, jpjm1 876 DO ji = 1, fs_jpim1 ! vector opt. 877 zwz(ji,jj) = ( ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 878 & - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) 879 END DO 880 END DO 881 CASE ( np_CRV ) !* Coriolis + relative vorticity 882 DO jj = 1, jpjm1 883 DO ji = 1, fs_jpim1 ! vector opt. 884 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 885 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 886 & * r1_e1e2f(ji,jj) ) 887 END DO 888 END DO 889 CASE ( np_CME ) !* Coriolis + metric 890 DO jj = 1, jpjm1 891 DO ji = 1, fs_jpim1 ! vector opt. 892 zwz(ji,jj) = ( ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 893 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) 894 END DO 895 END DO 896 CASE DEFAULT ! error 897 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 898 END SELECT 899 ! 900 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 901 DO jj = 1, jpjm1 902 DO ji = 1, fs_jpim1 ! vector opt. 903 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 904 END DO 905 END DO 906 ENDIF 907 ! 908 CALL lbc_lnk( zwz, 'F', 1. ) 909 ! 910 ! !== horizontal fluxes ==! 911 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 912 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 913 914 ! !== compute and add the vorticity term trend =! 915 jj = 2 916 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 917 DO ji = 2, jpi ! split in 2 parts due to vector opt. 918 ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) / ( e3u_n(ji , jj, jk) + e3v_n(ji, jj , jk ) ) 919 ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj , jk ) ) 920 ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) / ( e3u_n(ji , jj, jk) + e3v_n(ji, jj-1, jk ) ) 921 ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj-1, jk ) ) 922 END DO 923 DO jj = 3, jpj 924 DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3 925 ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) / ( e3u_n(ji , jj, jk) + e3v_n(ji, jj , jk ) ) 926 ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj , jk ) ) 927 ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) / ( e3u_n(ji , jj, jk) + e3v_n(ji, jj-1, jk ) ) 928 ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) / ( e3u_n(ji-1, jj, jk) + e3v_n(ji, jj-1, jk ) ) 929 END DO 930 END DO 931 DO jj = 2, jpjm1 932 DO ji = fs_2, fs_jpim1 ! vector opt. 933 zua = + r1_6 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 934 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 935 zva = - r1_6 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 936 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 937 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 938 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 939 END DO 940 END DO 941 ! ! =============== 942 END DO ! End of slab 943 ! ! =============== 944 END SUBROUTINE vor_eeUV 945 946 808 947 SUBROUTINE dyn_vor_init 809 948 !!--------------------------------------------------------------------- … … 817 956 !! 818 957 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 819 & ln_dynvor_ee n, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk958 & ln_dynvor_eeUV, ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 820 959 !!---------------------------------------------------------------------- 821 960 ! … … 840 979 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT 841 980 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT 981 WRITE(numout,*) ' energy conserving scheme (een using e3u and e3v) ln_dynvor_eeUV = ', ln_dynvor_eeUV 842 982 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 843 983 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f … … 873 1013 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 874 1014 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF 1015 IF( ln_dynvor_eeUV ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEUV ; ENDIF 875 1016 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 876 1017 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF … … 925 1066 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 926 1067 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 1068 CASE( np_EEUV ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3u and e3v) (EEUV)' 927 1069 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 928 1070 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)'
Note: See TracChangeset
for help on using the changeset viewer.
