Changeset 10874 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN
- Timestamp:
- 2019-04-15T15:57:37+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN
- Files:
-
- 9 edited
-
dynadv.F90 (modified) (3 diffs)
-
dynadv_cen2.F90 (modified) (9 diffs)
-
dynadv_ubs.F90 (modified) (10 diffs)
-
dynkeg.F90 (modified) (11 diffs)
-
dynldf.F90 (modified) (3 diffs)
-
dynldf_lap_blp.F90 (modified) (7 diffs)
-
dynvor.F90 (modified) (47 diffs)
-
dynzad.F90 (modified) (6 diffs)
-
dynzdf.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv.F90
r10789 r10874 53 53 CONTAINS 54 54 55 SUBROUTINE dyn_adv( kt , ktlev1, ktlev2, pu_rhs, pv_rhs)55 SUBROUTINE dyn_adv( kt ) 56 56 !!--------------------------------------------------------------------- 57 57 !! *** ROUTINE dyn_adv *** … … 67 67 !!---------------------------------------------------------------------- 68 68 INTEGER, INTENT( in ) :: kt ! ocean time-step index 69 INTEGER, INTENT( in ) :: ktlev1, ktlev2 ! time level indices for source terms70 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends71 69 !!---------------------------------------------------------------------- 72 70 ! … … 75 73 SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==! 76 74 CASE( np_VEC_c2 ) 77 CALL dyn_keg ( kt, ktlev2, nn_dynkeg, pu_rhs, pv_rhs) ! vector form : horizontal gradient of kinetic energy78 CALL dyn_zad ( kt , ktlev2, pu_rhs, pv_rhs )! vector form : vertical advection75 CALL dyn_keg ( kt, nn_dynkeg ) ! vector form : horizontal gradient of kinetic energy 76 CALL dyn_zad ( kt ) ! vector form : vertical advection 79 77 CASE( np_FLX_c2 ) 80 CALL dyn_adv_cen2( kt , ktlev2, pu_rhs, pv_rhs )! 2nd order centered scheme78 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme 81 79 CASE( np_FLX_ubs ) 82 CALL dyn_adv_ubs ( kt , ktlev1, ktlev2, pu_rhs, pv_rhs) ! 3rd order UBS scheme (UP3)80 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme (UP3) 83 81 END SELECT 84 82 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_cen2.F90
r10789 r10874 35 35 CONTAINS 36 36 37 SUBROUTINE dyn_adv_cen2( kt , ktlev, pu_rhs, pv_rhs)37 SUBROUTINE dyn_adv_cen2( kt ) 38 38 !!---------------------------------------------------------------------- 39 39 !! *** ROUTINE dyn_adv_cen2 *** … … 44 44 !! ** Method : Trend evaluated using now fields (centered in time) 45 45 !! 46 !! ** Action : ( pu_rhs,pv_rhs) updated with the now vorticity term trend46 !! ** Action : (ua,va) updated with the now vorticity term trend 47 47 !!---------------------------------------------------------------------- 48 INTEGER, INTENT( in ) :: kt ! ocean time-step index 49 INTEGER, INTENT( in ) :: ktlev ! time level index for source terms 50 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends 48 INTEGER, INTENT( in ) :: kt ! ocean time-step index 51 49 ! 52 50 INTEGER :: ji, jj, jk ! dummy loop indices … … 62 60 ! 63 61 IF( l_trddyn ) THEN ! trends: store the input trends 64 zfu_uw(:,:,:) = pu_rhs(:,:,:)65 zfv_vw(:,:,:) = pv_rhs(:,:,:)62 zfu_uw(:,:,:) = ua(:,:,:) 63 zfv_vw(:,:,:) = va(:,:,:) 66 64 ENDIF 67 65 ! … … 69 67 ! 70 68 DO jk = 1, jpkm1 ! horizontal transport 71 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u (:,:,jk,ktlev) * uu(:,:,jk,ktlev)72 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v (:,:,jk,ktlev) * vv(:,:,jk,ktlev)69 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 70 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 73 71 DO jj = 1, jpjm1 ! horizontal momentum fluxes (at T- and F-point) 74 72 DO ji = 1, fs_jpim1 ! vector opt. 75 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( u u(ji,jj,jk,ktlev) + uu(ji+1,jj ,jk,ktlev) )76 zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( u u(ji,jj,jk,ktlev) + uu(ji ,jj+1,jk,ktlev) )77 zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( v v(ji,jj,jk,ktlev) + vv(ji+1,jj ,jk,ktlev) )78 zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( v v(ji,jj,jk,ktlev) + vv(ji ,jj+1,jk,ktlev) )73 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) 74 zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji ,jj+1,jk) ) 75 zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji+1,jj ,jk) ) 76 zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) 79 77 END DO 80 78 END DO 81 79 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 82 80 DO ji = fs_2, fs_jpim1 ! vector opt. 83 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &84 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u (ji,jj,jk,ktlev)85 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &86 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v (ji,jj,jk,ktlev)81 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 82 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 83 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 84 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 87 85 END DO 88 86 END DO … … 90 88 ! 91 89 IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic 92 zfu_uw(:,:,:) = pu_rhs(:,:,:) - zfu_uw(:,:,:)93 zfv_vw(:,:,:) = pv_rhs(:,:,:) - zfv_vw(:,:,:)90 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 91 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 94 92 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 95 zfu_t(:,:,:) = pu_rhs(:,:,:)96 zfv_t(:,:,:) = pv_rhs(:,:,:)93 zfu_t(:,:,:) = ua(:,:,:) 94 zfv_t(:,:,:) = va(:,:,:) 97 95 ENDIF 98 96 ! … … 108 106 DO jj = 2, jpjm1 109 107 DO ji = fs_2, fs_jpim1 110 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * w w(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * uu(ji,jj,1,ktlev)111 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * w w(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * vv(ji,jj,1,ktlev)108 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 109 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 112 110 END DO 113 111 END DO … … 116 114 DO jj = 2, jpj ! 1/4 * Vertical transport 117 115 DO ji = 2, jpi 118 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * w w(ji,jj,jk)116 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 119 117 END DO 120 118 END DO 121 119 DO jj = 2, jpjm1 122 120 DO ji = fs_2, fs_jpim1 ! vector opt. 123 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( u u(ji,jj,jk,ktlev) + uu(ji,jj,jk-1,ktlev) )124 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( v v(ji,jj,jk,ktlev) + vv(ji,jj,jk-1,ktlev) )121 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 122 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 125 123 END DO 126 124 END DO … … 129 127 DO jj = 2, jpjm1 130 128 DO ji = fs_2, fs_jpim1 ! vector opt. 131 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)132 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)129 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 130 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 133 131 END DO 134 132 END DO … … 136 134 ! 137 135 IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic 138 zfu_t(:,:,:) = pu_rhs(:,:,:) - zfu_t(:,:,:)139 zfv_t(:,:,:) = pv_rhs(:,:,:) - zfv_t(:,:,:)136 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 137 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 140 138 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 141 139 ENDIF -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynadv_ubs.F90
r10802 r10874 41 41 CONTAINS 42 42 43 SUBROUTINE dyn_adv_ubs( kt , ktlev1, ktlev2, pu_rhs, pv_rhs)43 SUBROUTINE dyn_adv_ubs( kt ) 44 44 !!---------------------------------------------------------------------- 45 45 !! *** ROUTINE dyn_adv_ubs *** … … 64 64 !! gamma1=1/3 and gamma2=1/32. 65 65 !! 66 !! ** Action : - ( pu_rhs,pv_rhs) updated with the 3D advective momentum trends66 !! ** Action : - (ua,va) updated with the 3D advective momentum trends 67 67 !! 68 68 !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. 69 69 !!---------------------------------------------------------------------- 70 INTEGER, INTENT(in) :: kt ! ocean time-step index 71 INTEGER, INTENT(in) :: ktlev1, ktlev2 ! time level indices for source terms 72 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends 70 INTEGER, INTENT(in) :: kt ! ocean time-step index 73 71 ! 74 72 INTEGER :: ji, jj, jk ! dummy loop indices … … 97 95 ! 98 96 IF( l_trddyn ) THEN ! trends: store the input trends 99 zfu_uw(:,:,:) = pu_rhs(:,:,:)100 zfv_vw(:,:,:) = pv_rhs(:,:,:)97 zfu_uw(:,:,:) = ua(:,:,:) 98 zfv_vw(:,:,:) = va(:,:,:) 101 99 ENDIF 102 100 ! ! =========================== ! … … 104 102 ! ! =========================== ! 105 103 ! ! horizontal volume fluxes 106 zfu(:,:,jk) = e2u(:,:) * e3u (:,:,jk,ktlev2) * uu(:,:,jk,ktlev2)107 zfv(:,:,jk) = e1v(:,:) * e3v (:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)104 zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 105 zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 108 106 ! 109 107 DO jj = 2, jpjm1 ! laplacian 110 108 DO ji = fs_2, fs_jpim1 ! vector opt. 111 zlu_uu(ji,jj,jk,1) = ( u u (ji+1,jj ,jk,ktlev1) - 2.*uu (ji,jj,jk,ktlev1) + uu (ji-1,jj ,jk,ktlev1) ) * umask(ji,jj,jk)112 zlv_vv(ji,jj,jk,1) = ( v v (ji ,jj+1,jk,ktlev1) - 2.*vv (ji,jj,jk,ktlev1) + vv (ji ,jj-1,jk,ktlev1) ) * vmask(ji,jj,jk)113 zlu_uv(ji,jj,jk,1) = ( u u (ji ,jj+1,jk,ktlev1) - uu (ji ,jj ,jk,ktlev1) ) * fmask(ji ,jj ,jk) &114 & - ( u u (ji ,jj ,jk,ktlev1) - uu (ji ,jj-1,jk,ktlev1) ) * fmask(ji ,jj-1,jk)115 zlv_vu(ji,jj,jk,1) = ( v v (ji+1,jj ,jk,ktlev1) - vv (ji ,jj ,jk,ktlev1) ) * fmask(ji ,jj ,jk) &116 & - ( v v (ji ,jj ,jk,ktlev1) - vv (ji-1,jj ,jk,ktlev1) ) * fmask(ji-1,jj ,jk)109 zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj ,jk) - 2.*ub (ji,jj,jk) + ub (ji-1,jj ,jk) ) * umask(ji,jj,jk) 110 zlv_vv(ji,jj,jk,1) = ( vb (ji ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji ,jj-1,jk) ) * vmask(ji,jj,jk) 111 zlu_uv(ji,jj,jk,1) = ( ub (ji ,jj+1,jk) - ub (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 112 & - ( ub (ji ,jj ,jk) - ub (ji ,jj-1,jk) ) * fmask(ji ,jj-1,jk) 113 zlv_vu(ji,jj,jk,1) = ( vb (ji+1,jj ,jk) - vb (ji ,jj ,jk) ) * fmask(ji ,jj ,jk) & 114 & - ( vb (ji ,jj ,jk) - vb (ji-1,jj ,jk) ) * fmask(ji-1,jj ,jk) 117 115 ! 118 116 zlu_uu(ji,jj,jk,2) = ( zfu(ji+1,jj ,jk) - 2.*zfu(ji,jj,jk) + zfu(ji-1,jj ,jk) ) * umask(ji,jj,jk) … … 134 132 DO jk = 1, jpkm1 ! ====================== ! 135 133 ! ! horizontal volume fluxes 136 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u (:,:,jk,ktlev2) * uu(:,:,jk,ktlev2)137 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v (:,:,jk,ktlev2) * vv(:,:,jk,ktlev2)134 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 135 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 138 136 ! 139 137 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point 140 138 DO ji = 1, fs_jpim1 ! vector opt. 141 zui = ( u u(ji,jj,jk,ktlev2) + uu(ji+1,jj ,jk,ktlev2) )142 zvj = ( v v(ji,jj,jk,ktlev2) + vv(ji ,jj+1,jk,ktlev2) )139 zui = ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) 140 zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) 143 141 ! 144 142 IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) … … 166 164 ! 167 165 zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) & 168 & * ( u u(ji,jj,jk,ktlev2) + uu(ji ,jj+1,jk,ktlev2) - gamma1 * zl_u )166 & * ( un(ji,jj,jk) + un(ji ,jj+1,jk) - gamma1 * zl_u ) 169 167 zfu_f(ji ,jj ,jk) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,jk,2) + zlu_uv(ji ,jj+1,jk,2) ) ) & 170 & * ( v v(ji,jj,jk,ktlev2) + vv(ji+1,jj ,jk,ktlev2) - gamma1 * zl_v )168 & * ( vn(ji,jj,jk) + vn(ji+1,jj ,jk) - gamma1 * zl_v ) 171 169 END DO 172 170 END DO 173 171 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 174 172 DO ji = fs_2, fs_jpim1 ! vector opt. 175 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &176 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u (ji,jj,jk,ktlev2)177 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &178 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v (ji,jj,jk,ktlev2)173 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 174 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 175 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 176 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 179 177 END DO 180 178 END DO 181 179 END DO 182 180 IF( l_trddyn ) THEN ! trends: send trends to trddyn for diagnostic 183 zfu_uw(:,:,:) = pu_rhs(:,:,:) - zfu_uw(:,:,:)184 zfv_vw(:,:,:) = pv_rhs(:,:,:) - zfv_vw(:,:,:)181 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 182 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 185 183 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 186 zfu_t(:,:,:) = pu_rhs(:,:,:)187 zfv_t(:,:,:) = pv_rhs(:,:,:)184 zfu_t(:,:,:) = ua(:,:,:) 185 zfv_t(:,:,:) = va(:,:,:) 188 186 ENDIF 189 187 ! ! ==================== ! … … 201 199 DO jj = 2, jpjm1 202 200 DO ji = fs_2, fs_jpim1 203 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * w w(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * uu(ji,jj,1,ktlev2)204 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * w w(ji,jj,1) + e1e2t(ji,jj+1) * ww(ji,jj+1,1) ) * vv(ji,jj,1,ktlev2)201 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 202 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 205 203 END DO 206 204 END DO … … 209 207 DO jj = 2, jpj 210 208 DO ji = 2, jpi 211 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * w w(ji,jj,jk)209 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 212 210 END DO 213 211 END DO 214 212 DO jj = 2, jpjm1 215 213 DO ji = fs_2, fs_jpim1 ! vector opt. 216 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( u u(ji,jj,jk,ktlev2) + uu(ji,jj,jk-1,ktlev2) )217 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( v v(ji,jj,jk,ktlev2) + vv(ji,jj,jk-1,ktlev2) )214 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 215 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 218 216 END DO 219 217 END DO … … 222 220 DO jj = 2, jpjm1 223 221 DO ji = fs_2, fs_jpim1 ! vector opt. 224 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)225 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev2)222 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 223 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 226 224 END DO 227 225 END DO … … 229 227 ! 230 228 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 231 zfu_t(:,:,:) = pu_rhs(:,:,:) - zfu_t(:,:,:)232 zfv_t(:,:,:) = pv_rhs(:,:,:) - zfv_t(:,:,:)229 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 230 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 233 231 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 234 232 ENDIF -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynkeg.F90
r10789 r10874 44 44 CONTAINS 45 45 46 SUBROUTINE dyn_keg( kt, k tlev, kscheme, pu_rhs, pv_rhs)46 SUBROUTINE dyn_keg( kt, kscheme ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE dyn_keg *** … … 54 54 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 55 55 !! conserve kinetic energy. Compute the now horizontal kinetic energy 56 !! zhke = 1/2 [ mi-1( u u(:,:,:,ktlev)^2 ) + mj-1( vv(:,:,:,ktlev)^2 ) ]56 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 57 57 !! * kscheme = nkeg_HW : Hollingsworth correction following 58 58 !! Arakawa (2001). The now horizontal kinetic energy is given by: 59 !! zhke = 1/6 [ mi-1( 2 * u u(:,:,:,ktlev)^2 + ((uu(j+1,ktlev)+uu(j-1,ktlev))/2)^2 )60 !! + mj-1( 2 * v v(:,:,:,ktlev)^2 + ((vv(i+1,ktlev)+vv(i-1,ktlev))/2)^2 ) ]59 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((un(j+1)+un(j-1))/2)^2 ) 60 !! + mj-1( 2 * vn^2 + ((vn(i+1)+vn(i-1))/2)^2 ) ] 61 61 !! 62 62 !! Take its horizontal gradient and add it to the general momentum 63 !! trend ( pu_rhs,pv_rhs).64 !! pu_rhs = pu_rhs- 1/e1u di[ zhke ]65 !! pv_rhs = pv_rhs- 1/e2v dj[ zhke ]63 !! trend (ua,va). 64 !! ua = ua - 1/e1u di[ zhke ] 65 !! va = va - 1/e2v dj[ zhke ] 66 66 !! 67 !! ** Action : - Update the ( pu_rhs, pv_rhs) with the hor. ke gradient trend67 !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 68 68 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 69 69 !! … … 72 72 !!---------------------------------------------------------------------- 73 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 INTEGER, INTENT( in ) :: ktlev ! time level index for source terms75 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 76 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends77 75 ! 78 76 INTEGER :: ji, jj, jk, jb ! dummy loop indices … … 94 92 IF( l_trddyn ) THEN ! Save the input trends 95 93 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 96 ztrdu(:,:,:) = pu_rhs(:,:,:)97 ztrdv(:,:,:) = pv_rhs(:,:,:)94 ztrdu(:,:,:) = ua(:,:,:) 95 ztrdv(:,:,:) = va(:,:,:) 98 96 ENDIF 99 97 … … 111 109 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 112 110 ifu = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 113 u u(ii-ifu,ij,jk,ktlev) = uu(ii,ij,jk,ktlev) * umask(ii,ij,jk)111 un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 114 112 END DO 115 113 END DO … … 121 119 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 122 120 ifv = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 123 v v(ii,ij-ifv,jk,ktlev) = vv(ii,ij,jk,ktlev) * vmask(ii,ij,jk)121 vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 124 122 END DO 125 123 END DO … … 134 132 DO jj = 2, jpj 135 133 DO ji = fs_2, jpi ! vector opt. 136 zu = u u(ji-1,jj ,jk,ktlev) * uu(ji-1,jj ,jk,ktlev) &137 & + u u(ji ,jj ,jk,ktlev) * uu(ji ,jj ,jk,ktlev)138 zv = v v(ji ,jj-1,jk,ktlev) * vv(ji ,jj-1,jk,ktlev) &139 & + v v(ji ,jj ,jk,ktlev) * vv(ji ,jj ,jk,ktlev)134 zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 135 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) 136 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 137 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) 140 138 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 141 139 END DO … … 147 145 DO jj = 2, jpjm1 148 146 DO ji = fs_2, jpim1 ! vector opt. 149 zu = 8._wp * ( u u(ji-1,jj ,jk,ktlev) * uu(ji-1,jj ,jk,ktlev) &150 & + u u(ji ,jj ,jk,ktlev) * uu(ji ,jj ,jk,ktlev) ) &151 & + ( u u(ji-1,jj-1,jk,ktlev) + uu(ji-1,jj+1,jk,ktlev) ) * ( uu(ji-1,jj-1,jk,ktlev) + uu(ji-1,jj+1,jk,ktlev) ) &152 & + ( u u(ji ,jj-1,jk,ktlev) + uu(ji ,jj+1,jk,ktlev) ) * ( uu(ji ,jj-1,jk,ktlev) + uu(ji ,jj+1,jk,ktlev) )147 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 148 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) & 149 & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) & 150 & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) 153 151 ! 154 zv = 8._wp * ( v v(ji ,jj-1,jk,ktlev) * vv(ji ,jj-1,jk,ktlev) &155 & + v v(ji ,jj ,jk,ktlev) * vv(ji ,jj ,jk,ktlev) ) &156 & + ( v v(ji-1,jj-1,jk,ktlev) + vv(ji+1,jj-1,jk,ktlev) ) * ( vv(ji-1,jj-1,jk,ktlev) + vv(ji+1,jj-1,jk,ktlev) ) &157 & + ( v v(ji-1,jj ,jk,ktlev) + vv(ji+1,jj ,jk,ktlev) ) * ( vv(ji-1,jj ,jk,ktlev) + vv(ji+1,jj ,jk,ktlev) )152 zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 153 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) & 154 & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) & 155 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) 158 156 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 159 157 END DO … … 166 164 IF (ln_bdy) THEN 167 165 ! restore velocity masks at points outside boundary 168 u u(:,:,:,ktlev) = uu(:,:,:,ktlev) * umask(:,:,:)169 v v(:,:,:,ktlev) = vv(:,:,:,ktlev) * vmask(:,:,:)166 un(:,:,:) = un(:,:,:) * umask(:,:,:) 167 vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 170 168 ENDIF 171 169 … … 174 172 DO jj = 2, jpjm1 175 173 DO ji = fs_2, fs_jpim1 ! vector opt. 176 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj)177 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj)174 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 175 va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 178 176 END DO 179 177 END DO … … 181 179 ! 182 180 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 183 ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)184 ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)181 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 182 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 185 183 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 186 184 DEALLOCATE( ztrdu , ztrdv ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf.F90
r10806 r10874 43 43 CONTAINS 44 44 45 SUBROUTINE dyn_ldf( kt , ktlev1, ktlev2, pu_rhs, pv_rhs)45 SUBROUTINE dyn_ldf( kt ) 46 46 !!---------------------------------------------------------------------- 47 47 !! *** ROUTINE dyn_ldf *** … … 49 49 !! ** Purpose : compute the lateral ocean dynamics physics. 50 50 !!---------------------------------------------------------------------- 51 INTEGER, INTENT(in) :: kt ! ocean time-step index 52 INTEGER, INTENT(in) :: ktlev1, ktlev2 ! time level index for source terms 53 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends 51 INTEGER, INTENT(in) :: kt ! ocean time-step index 54 52 ! 55 53 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv … … 60 58 IF( l_trddyn ) THEN ! temporary save of momentum trends 61 59 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 62 ztrdu(:,:,:) = pu_rhs(:,:,:)63 ztrdv(:,:,:) = pv_rhs(:,:,:)60 ztrdu(:,:,:) = ua(:,:,:) 61 ztrdv(:,:,:) = va(:,:,:) 64 62 ENDIF 65 63 66 64 SELECT CASE ( nldf_dyn ) ! compute lateral mixing trend and add it to the general trend 67 65 ! 68 CASE ( np_lap ) ; CALL dyn_ldf_lap( kt, ktlev1, ktlev2, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs, 1 ) ! iso-level laplacian66 CASE ( np_lap ) ; CALL dyn_ldf_lap( kt, ub, vb, ua, va, 1 ) ! iso-level laplacian 69 67 CASE ( np_lap_i ) ; CALL dyn_ldf_iso( kt ) ! rotated laplacian 70 CASE ( np_blp ) ; CALL dyn_ldf_blp( kt, ktlev1, ktlev2, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs) ! iso-level bi-laplacian68 CASE ( np_blp ) ; CALL dyn_ldf_blp( kt, ub, vb, ua, va ) ! iso-level bi-laplacian 71 69 ! 72 70 END SELECT 73 71 74 72 IF( l_trddyn ) THEN ! save the horizontal diffusive trends for further diagnostics 75 ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)76 ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)73 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 74 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 77 75 CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 78 76 DEALLOCATE ( ztrdu , ztrdv ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90
r10806 r10874 35 35 CONTAINS 36 36 37 SUBROUTINE dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs, kpass )37 SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass ) 38 38 !!---------------------------------------------------------------------- 39 39 !! *** ROUTINE dyn_ldf_lap *** … … 45 45 !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 46 46 !! 47 !! ** Action : - pu _rhs, pva_rhs increased by the harmonic operator applied on pu, pv.47 !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb. 48 48 !!---------------------------------------------------------------------- 49 INTEGER , INTENT(in ) :: kt ! ocean time-step index 50 INTEGER , INTENT(in ) :: ktlev1, ktlev2 ! time level index for scale factors 49 INTEGER , INTENT(in ) :: kt ! ocean time-step index 51 50 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 52 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv! before velocity [m/s]53 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu _rhs, pva_rhs! velocity trend [m/s2]51 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity [m/s] 52 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! velocity trend [m/s2] 54 53 ! 55 54 INTEGER :: ji, jj, jk ! dummy loop indices … … 77 76 !!gm open question here : e3f at before or now ? probably now... 78 77 !!gm note that ahmf has already been multiplied by fmask 79 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f (ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) &80 & * ( e2v(ji ,jj-1) * pv (ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) &81 & - e1u(ji-1,jj ) * pu (ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) )78 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 79 & * ( e2v(ji ,jj-1) * pvb(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk) & 80 & - e1u(ji-1,jj ) * pub(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk) ) 82 81 ! ! ahm * div (computed from 2 to jpi/jpj) 83 82 !!gm note that ahmt has already been multiplied by tmask 84 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t (ji,jj,jk,ktlev1) &85 & * ( e2u(ji,jj)*e3u (ji,jj,jk,ktlev1) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,ktlev1) * pu(ji-1,jj,jk) &86 & + e1v(ji,jj)*e3v (ji,jj,jk,ktlev1) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,ktlev1) * pv(ji,jj-1,jk) )83 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) & 84 & * ( e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk) & 85 & + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk) ) 87 86 END DO 88 87 END DO … … 90 89 DO jj = 2, jpjm1 ! - curl( curl) + grad( div ) 91 90 DO ji = fs_2, fs_jpim1 ! vector opt. 92 pu _rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( &93 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u (ji,jj,jk,ktlev2) &91 pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * ( & 92 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) & 94 93 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 95 94 ! 96 pva _rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zsign * ( &97 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v (ji,jj,jk,ktlev2) &95 pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * ( & 96 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) & 98 97 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 99 98 END DO … … 106 105 107 106 108 SUBROUTINE dyn_ldf_blp( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs)107 SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva ) 109 108 !!---------------------------------------------------------------------- 110 109 !! *** ROUTINE dyn_ldf_blp *** … … 117 116 !! It is computed by two successive calls to dyn_ldf_lap routine 118 117 !! 119 !! ** Action : pt _rhsupdated with the before rotated bilaplacian diffusion118 !! ** Action : pta updated with the before rotated bilaplacian diffusion 120 119 !!---------------------------------------------------------------------- 121 120 INTEGER , INTENT(in ) :: kt ! ocean time-step index 122 INTEGER , INTENT(in ) :: ktlev1, ktlev2 ! time level index for scale factors 123 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv ! before velocity fields 124 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pva_rhs ! momentum trend 121 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pub, pvb ! before velocity fields 122 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! momentum trend 125 123 ! 126 124 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulap, zvlap ! laplacian at u- and v-point … … 136 134 zvlap(:,:,:) = 0._wp 137 135 ! 138 CALL dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, zulap, zvlap, 1 ) ! rotated laplacian applied to pt(output in zlap)136 CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 ) ! rotated laplacian applied to ptb (output in zlap) 139 137 ! 140 138 CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. ) ! Lateral boundary conditions 141 139 ! 142 CALL dyn_ldf_lap( kt, ktlev1, ktlev2, zulap, zvlap, pu_rhs, pva_rhs, 2 ) ! rotated laplacian applied to zlap (output in pt_rhs)140 CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 ) ! rotated laplacian applied to zlap (output in pta) 143 141 ! 144 142 END SUBROUTINE dyn_ldf_blp -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynvor.F90
r10806 r10874 96 96 CONTAINS 97 97 98 SUBROUTINE dyn_vor( kt , ktlev, pu_rhs, pv_rhs)98 SUBROUTINE dyn_vor( kt ) 99 99 !!---------------------------------------------------------------------- 100 100 !! 101 101 !! ** Purpose : compute the lateral ocean tracer physics. 102 102 !! 103 !! ** Action : - Update ( pu_rhs,pv_rhs) with the now vorticity term trend103 !! ** Action : - Update (ua,va) with the now vorticity term trend 104 104 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative 105 105 !! and planetary vorticity trends) and send them to trd_dyn 106 106 !! for futher diagnostics (l_trddyn=T) 107 107 !!---------------------------------------------------------------------- 108 INTEGER, INTENT( in ) :: kt ! ocean time-step index 109 INTEGER, INTENT( in ) :: ktlev ! time level index for source terms 110 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends 108 INTEGER, INTENT( in ) :: kt ! ocean time-step index 111 109 ! 112 110 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv … … 119 117 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 120 118 ! 121 ztrdu(:,:,:) = pu_rhs(:,:,:) !* planetary vorticity trend (including Stokes-Coriolis force)122 ztrdv(:,:,:) = pv_rhs(:,:,:)119 ztrdu(:,:,:) = ua(:,:,:) !* planetary vorticity trend (including Stokes-Coriolis force) 120 ztrdv(:,:,:) = va(:,:,:) 123 121 SELECT CASE( nvor_scheme ) 124 CASE( np_ENS ) ; CALL vor_ens( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! enstrophy conserving scheme125 IF( ln_stcor ) CALL vor_ens( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend126 CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy conserving scheme127 IF( ln_stcor ) CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend128 CASE( np_ENT ) ; CALL vor_enT( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy conserving scheme (T-pts)129 IF( ln_stcor ) CALL vor_enT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend130 CASE( np_EET ) ; CALL vor_eeT( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy conserving scheme (een with e3t)131 IF( ln_stcor ) CALL vor_eeT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend132 CASE( np_EEN ) ; CALL vor_een( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy & enstrophy scheme133 IF( ln_stcor ) CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend122 CASE( np_ENS ) ; CALL vor_ens( kt, ncor, un , vn , ua, va ) ! enstrophy conserving scheme 123 IF( ln_stcor ) CALL vor_ens( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 124 CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, ncor, un , vn , ua, va ) ! energy conserving scheme 125 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 126 CASE( np_ENT ) ; CALL vor_enT( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (T-pts) 127 IF( ln_stcor ) CALL vor_enT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 128 CASE( np_EET ) ; CALL vor_eeT( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 129 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 130 CASE( np_EEN ) ; CALL vor_een( kt, ncor, un , vn , ua, va ) ! energy & enstrophy scheme 131 IF( ln_stcor ) CALL vor_een( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 134 132 END SELECT 135 ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)136 ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)133 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 134 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 137 135 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 138 136 ! 139 137 IF( n_dynadv /= np_LIN_dyn ) THEN !* relative vorticity or metric trend (only in non-linear case) 140 ztrdu(:,:,:) = pu_rhs(:,:,:)141 ztrdv(:,:,:) = pv_rhs(:,:,:)138 ztrdu(:,:,:) = ua(:,:,:) 139 ztrdv(:,:,:) = va(:,:,:) 142 140 SELECT CASE( nvor_scheme ) 143 CASE( np_ENT ) ; CALL vor_enT( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy conserving scheme (T-pts)144 CASE( np_EET ) ; CALL vor_eeT( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy conserving scheme (een with e3t)145 CASE( np_ENE ) ; CALL vor_ene( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy conserving scheme146 CASE( np_ENS, np_MIX ) ; CALL vor_ens( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! enstrophy conserving scheme147 CASE( np_EEN ) ; CALL vor_een( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! energy & enstrophy scheme141 CASE( np_ENT ) ; CALL vor_enT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (T-pts) 142 CASE( np_EET ) ; CALL vor_eeT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 143 CASE( np_ENE ) ; CALL vor_ene( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme 144 CASE( np_ENS, np_MIX ) ; CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! enstrophy conserving scheme 145 CASE( np_EEN ) ; CALL vor_een( kt, nrvm, un , vn , ua, va ) ! energy & enstrophy scheme 148 146 END SELECT 149 ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)150 ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)147 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 148 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 151 149 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 152 150 ENDIF … … 158 156 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 159 157 CASE( np_ENT ) !* energy conserving scheme (T-pts) 160 CALL vor_enT( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! total vorticity trend161 IF( ln_stcor ) CALL vor_enT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend158 CALL vor_enT( kt, ntot, un , vn , ua, va ) ! total vorticity trend 159 IF( ln_stcor ) CALL vor_enT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 162 160 CASE( np_EET ) !* energy conserving scheme (een scheme using e3t) 163 CALL vor_eeT( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! total vorticity trend164 IF( ln_stcor ) CALL vor_eeT( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend161 CALL vor_eeT( kt, ntot, un , vn , ua, va ) ! total vorticity trend 162 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 165 163 CASE( np_ENE ) !* energy conserving scheme 166 CALL vor_ene( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! total vorticity trend167 IF( ln_stcor ) CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend164 CALL vor_ene( kt, ntot, un , vn , ua, va ) ! total vorticity trend 165 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 168 166 CASE( np_ENS ) !* enstrophy conserving scheme 169 CALL vor_ens( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! total vorticity trend170 IF( ln_stcor ) CALL vor_ens( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend167 CALL vor_ens( kt, ntot, un , vn , ua, va ) ! total vorticity trend 168 IF( ln_stcor ) CALL vor_ens( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 171 169 CASE( np_MIX ) !* mixed ene-ens scheme 172 CALL vor_ens( kt, ktlev, nrvm, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! relative vorticity or metric trend (ens)173 CALL vor_ene( kt, ktlev, ncor, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! planetary vorticity trend (ene)174 IF( ln_stcor ) CALL vor_ene( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend170 CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! relative vorticity or metric trend (ens) 171 CALL vor_ene( kt, ncor, un , vn , ua, va ) ! planetary vorticity trend (ene) 172 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 175 173 CASE( np_EEN ) !* energy and enstrophy conserving scheme 176 CALL vor_een( kt, ktlev, ntot, uu(:,:,:,ktlev) , vv(:,:,:,ktlev) , pu_rhs, pv_rhs) ! total vorticity trend177 IF( ln_stcor ) CALL vor_een( kt, ktlev, ncor, usd, vsd, pu_rhs, pv_rhs) ! add the Stokes-Coriolis trend174 CALL vor_een( kt, ntot, un , vn , ua, va ) ! total vorticity trend 175 IF( ln_stcor ) CALL vor_een( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 178 176 END SELECT 179 177 ! … … 189 187 190 188 191 SUBROUTINE vor_enT( kt, k tlev, kvor, pu, pv, pu_rhs, pv_rhs )189 SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs ) 192 190 !!---------------------------------------------------------------------- 193 191 !! *** ROUTINE vor_enT *** … … 205 203 !! where rvor is the relative vorticity at f-point 206 204 !! 207 !! ** Action : - Update (u _rhs,v_rhs) with the now vorticity term trend205 !! ** Action : - Update (ua,va) with the now vorticity term trend 208 206 !!---------------------------------------------------------------------- 209 207 INTEGER , INTENT(in ) :: kt ! ocean time-step index 210 INTEGER , INTENT( in ) :: ktlev ! time level index for source terms211 208 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 212 209 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities … … 273 270 SELECT CASE( kvor ) !== volume weighted vorticity considered ==! 274 271 CASE ( np_COR ) !* Coriolis (planetary vorticity) 275 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t (:,:,jk,ktlev)272 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 276 273 CASE ( np_RVO ) !* relative vorticity 277 274 DO jj = 2, jpj 278 275 DO ji = 2, jpi ! vector opt. 279 276 zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 280 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t (ji,jj,jk,ktlev)277 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 281 278 END DO 282 279 END DO … … 285 282 DO ji = 2, jpi 286 283 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 287 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t (ji,jj,jk,ktlev)284 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t_n(ji,jj,jk) 288 285 END DO 289 286 END DO … … 292 289 DO ji = 2, jpi ! vector opt. 293 290 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 294 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) * e1e2t(ji,jj)*e3t (ji,jj,jk,ktlev)291 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 295 292 END DO 296 293 END DO … … 300 297 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 301 298 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 302 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t (ji,jj,jk,ktlev)299 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t_n(ji,jj,jk) 303 300 END DO 304 301 END DO … … 310 307 DO jj = 2, jpjm1 311 308 DO ji = 2, jpim1 ! vector opt. 312 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u (ji,jj,jk,ktlev) &309 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) & 313 310 & * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) & 314 311 & + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) ) 315 312 ! 316 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v (ji,jj,jk,ktlev) &313 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) & 317 314 & * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) & 318 315 & + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) ) … … 325 322 326 323 327 SUBROUTINE vor_ene( kt, k tlev, kvor, pu, pv, pu_rhs, pva_rhs)324 SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 328 325 !!---------------------------------------------------------------------- 329 326 !! *** ROUTINE vor_ene *** … … 337 334 !! The general trend of momentum is increased due to the vorticity 338 335 !! term which is given by: 339 !! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v v v(:,:,:,ktlev)) ]340 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u u u(:,:,:,ktlev)) ]336 !! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v vn) ] 337 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u un) ] 341 338 !! where rvor is the relative vorticity 342 339 !! 343 !! ** Action : - Update (u _rhs,v_rhs) with the now vorticity term trend340 !! ** Action : - Update (ua,va) with the now vorticity term trend 344 341 !! 345 342 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 346 343 !!---------------------------------------------------------------------- 347 344 INTEGER , INTENT(in ) :: kt ! ocean time-step index 348 INTEGER , INTENT( in ) :: ktlev ! time level index for source terms349 345 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 350 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu , pv! now velocities351 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu _rhs, pva_rhs! total v-trend346 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun, pvn ! now velocities 347 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! total v-trend 352 348 ! 353 349 INTEGER :: ji, jj, jk ! dummy loop indices … … 372 368 DO jj = 1, jpjm1 373 369 DO ji = 1, fs_jpim1 ! vector opt. 374 zwz(ji,jj) = ( e2v(ji+1,jj ) * pv (ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &375 & - e1u(ji ,jj+1) * pu (ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)370 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 371 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 376 372 END DO 377 373 END DO … … 379 375 DO jj = 1, jpjm1 380 376 DO ji = 1, fs_jpim1 ! vector opt. 381 zwz(ji,jj) = ( pv (ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &382 & - ( pu (ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)377 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 378 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 383 379 END DO 384 380 END DO … … 386 382 DO jj = 1, jpjm1 387 383 DO ji = 1, fs_jpim1 ! vector opt. 388 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pv (ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &389 & - e1u(ji,jj+1) * pu (ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)384 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 385 & - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 390 386 END DO 391 387 END DO … … 393 389 DO jj = 1, jpjm1 394 390 DO ji = 1, fs_jpim1 ! vector opt. 395 zwz(ji,jj) = ff_f(ji,jj) + ( pv (ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &396 & - ( pu (ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)391 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 392 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 397 393 END DO 398 394 END DO … … 410 406 411 407 IF( ln_sco ) THEN 412 zwz(:,:) = zwz(:,:) / e3f (:,:,jk)413 zwx(:,:) = e2u(:,:) * e3u (:,:,jk,ktlev) * pu(:,:,jk)414 zwy(:,:) = e1v(:,:) * e3v (:,:,jk,ktlev) * pv(:,:,jk)408 zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 409 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 410 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 415 411 ELSE 416 zwx(:,:) = e2u(:,:) * pu (:,:,jk)417 zwy(:,:) = e1v(:,:) * pv (:,:,jk)412 zwx(:,:) = e2u(:,:) * pun(:,:,jk) 413 zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 418 414 ENDIF 419 415 ! !== compute and add the vorticity term trend =! … … 424 420 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 425 421 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 426 pu _rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )427 pva _rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )422 pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 423 pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 428 424 END DO 429 425 END DO … … 434 430 435 431 436 SUBROUTINE vor_ens( kt, k tlev, kvor, pu, pv, pu_rhs, pva_rhs)432 SUBROUTINE vor_ens( kt, kvor, pun, pvn, pua, pva ) 437 433 !!---------------------------------------------------------------------- 438 434 !! *** ROUTINE vor_ens *** … … 445 441 !! potential enstrophy of a horizontally non-divergent flow. the 446 442 !! trend of the vorticity term is given by: 447 !! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v v v(:,:,:,ktlev)) ]448 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u u u(:,:,:,ktlev)) ]449 !! Add this trend to the general momentum trend (u _rhs,v_rhs):450 !! (u _rhs,v_rhs) = (u_rhs,v_rhs) + ( voru , vorv )451 !! 452 !! ** Action : - Update (u _rhs,v_rhs) arrays with the now vorticity term trend443 !! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 444 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 445 !! Add this trend to the general momentum trend (ua,va): 446 !! (ua,va) = (ua,va) + ( voru , vorv ) 447 !! 448 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 453 449 !! 454 450 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 455 451 !!---------------------------------------------------------------------- 456 452 INTEGER , INTENT(in ) :: kt ! ocean time-step index 457 INTEGER , INTENT( in ) :: ktlev ! time level index for source terms458 453 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 459 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu , pv! now velocities460 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu _rhs, pva_rhs! total v-trend454 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun, pvn ! now velocities 455 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! total v-trend 461 456 ! 462 457 INTEGER :: ji, jj, jk ! dummy loop indices … … 480 475 DO jj = 1, jpjm1 481 476 DO ji = 1, fs_jpim1 ! vector opt. 482 zwz(ji,jj) = ( e2v(ji+1,jj ) * pv (ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &483 & - e1u(ji ,jj+1) * pu (ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)477 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 478 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 484 479 END DO 485 480 END DO … … 487 482 DO jj = 1, jpjm1 488 483 DO ji = 1, fs_jpim1 ! vector opt. 489 zwz(ji,jj) = ( pv (ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &490 & - ( pu (ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)484 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 485 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 491 486 END DO 492 487 END DO … … 494 489 DO jj = 1, jpjm1 495 490 DO ji = 1, fs_jpim1 ! vector opt. 496 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv (ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &497 & - e1u(ji ,jj+1) * pu (ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)491 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 492 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 498 493 END DO 499 494 END DO … … 501 496 DO jj = 1, jpjm1 502 497 DO ji = 1, fs_jpim1 ! vector opt. 503 zwz(ji,jj) = ff_f(ji,jj) + ( pv (ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &504 & - ( pu (ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)498 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 499 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 505 500 END DO 506 501 END DO … … 518 513 ! 519 514 IF( ln_sco ) THEN !== horizontal fluxes ==! 520 zwz(:,:) = zwz(:,:) / e3f (:,:,jk)521 zwx(:,:) = e2u(:,:) * e3u (:,:,jk,ktlev) * pu(:,:,jk)522 zwy(:,:) = e1v(:,:) * e3v (:,:,jk,ktlev) * pv(:,:,jk)515 zwz(:,:) = zwz(:,:) / e3f_n(:,:,jk) 516 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 517 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 523 518 ELSE 524 zwx(:,:) = e2u(:,:) * pu (:,:,jk)525 zwy(:,:) = e1v(:,:) * pv (:,:,jk)519 zwx(:,:) = e2u(:,:) * pun(:,:,jk) 520 zwy(:,:) = e1v(:,:) * pvn(:,:,jk) 526 521 ENDIF 527 522 ! !== compute and add the vorticity term trend =! … … 532 527 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 533 528 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 534 pu _rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) )535 pva _rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) )529 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 530 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 536 531 END DO 537 532 END DO … … 542 537 543 538 544 SUBROUTINE vor_een( kt, k tlev, kvor, pu, pv, pu_rhs, pva_rhs)539 SUBROUTINE vor_een( kt, kvor, pun, pvn, pua, pva ) 545 540 !!---------------------------------------------------------------------- 546 541 !! *** ROUTINE vor_een *** … … 553 548 !! both the horizontal kinetic energy and the potential enstrophy 554 549 !! when horizontal divergence is zero (see the NEMO documentation) 555 !! Add this trend to the general momentum trend (u _rhs,v_rhs).556 !! 557 !! ** Action : - Update (u _rhs,v_rhs) with the now vorticity term trend550 !! Add this trend to the general momentum trend (ua,va). 551 !! 552 !! ** Action : - Update (ua,va) with the now vorticity term trend 558 553 !! 559 554 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 560 555 !!---------------------------------------------------------------------- 561 556 INTEGER , INTENT(in ) :: kt ! ocean time-step index 562 INTEGER , INTENT( in ) :: ktlev ! time level index for source terms563 557 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 564 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu , pv! now velocities565 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu _rhs, pva_rhs! total v-trend558 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun, pvn ! now velocities 559 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! total v-trend 566 560 ! 567 561 INTEGER :: ji, jj, jk ! dummy loop indices … … 588 582 DO jj = 1, jpjm1 589 583 DO ji = 1, fs_jpim1 ! vector opt. 590 ze3f = ( e3t (ji,jj+1,jk,ktlev)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,ktlev)*tmask(ji+1,jj+1,jk) &591 & + e3t (ji,jj ,jk,ktlev)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,ktlev)*tmask(ji+1,jj ,jk) )584 ze3f = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 585 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 592 586 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f 593 587 ELSE ; z1_e3f(ji,jj) = 0._wp … … 598 592 DO jj = 1, jpjm1 599 593 DO ji = 1, fs_jpim1 ! vector opt. 600 ze3f = ( e3t (ji,jj+1,jk,ktlev)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,ktlev)*tmask(ji+1,jj+1,jk) &601 & + e3t (ji,jj ,jk,ktlev)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,ktlev)*tmask(ji+1,jj ,jk) )594 ze3f = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 595 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 602 596 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 603 597 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) … … 619 613 DO jj = 1, jpjm1 620 614 DO ji = 1, fs_jpim1 ! vector opt. 621 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv (ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &622 & - e1u(ji ,jj+1) * pu (ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)615 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 616 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 623 617 END DO 624 618 END DO … … 626 620 DO jj = 1, jpjm1 627 621 DO ji = 1, fs_jpim1 ! vector opt. 628 zwz(ji,jj,jk) = ( ( pv (ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &629 & - ( pu (ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)622 zwz(ji,jj,jk) = ( ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 623 & - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 630 624 END DO 631 625 END DO … … 633 627 DO jj = 1, jpjm1 634 628 DO ji = 1, fs_jpim1 ! vector opt. 635 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv (ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &636 & - e1u(ji ,jj+1) * pu (ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &629 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 630 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 637 631 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 638 632 END DO … … 641 635 DO jj = 1, jpjm1 642 636 DO ji = 1, fs_jpim1 ! vector opt. 643 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv (ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &644 & - ( pu (ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)637 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 638 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 645 639 END DO 646 640 END DO … … 663 657 ! 664 658 ! !== horizontal fluxes ==! 665 zwx(:,:) = e2u(:,:) * e3u (:,:,jk,ktlev) * pu(:,:,jk)666 zwy(:,:) = e1v(:,:) * e3v (:,:,jk,ktlev) * pv(:,:,jk)659 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 660 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 667 661 668 662 ! !== compute and add the vorticity term trend =! … … 689 683 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 690 684 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 691 pu _rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua692 pva _rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zva685 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 686 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 693 687 END DO 694 688 END DO … … 700 694 701 695 702 SUBROUTINE vor_eeT( kt, k tlev, kvor, pu, pv, pu_rhs, pva_rhs)696 SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva ) 703 697 !!---------------------------------------------------------------------- 704 698 !! *** ROUTINE vor_eeT *** … … 711 705 !! a modified version of Arakawa and Lamb (1980) scheme (see vor_een). 712 706 !! The change consists in 713 !! Add this trend to the general momentum trend (u _rhs,v_rhs).714 !! 715 !! ** Action : - Update (u _rhs,v_rhs) with the now vorticity term trend707 !! Add this trend to the general momentum trend (ua,va). 708 !! 709 !! ** Action : - Update (ua,va) with the now vorticity term trend 716 710 !! 717 711 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 718 712 !!---------------------------------------------------------------------- 719 713 INTEGER , INTENT(in ) :: kt ! ocean time-step index 720 INTEGER , INTENT( in ) :: ktlev ! time level index for source terms721 714 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 722 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu , pv! now velocities723 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu _rhs, pva_rhs! total v-trend715 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun, pvn ! now velocities 716 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! total v-trend 724 717 ! 725 718 INTEGER :: ji, jj, jk ! dummy loop indices … … 753 746 DO jj = 1, jpjm1 754 747 DO ji = 1, fs_jpim1 ! vector opt. 755 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv (ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &756 & - e1u(ji ,jj+1) * pu (ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &748 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 749 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 757 750 & * r1_e1e2f(ji,jj) 758 751 END DO … … 761 754 DO jj = 1, jpjm1 762 755 DO ji = 1, fs_jpim1 ! vector opt. 763 zwz(ji,jj,jk) = ( pv (ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &764 & - ( pu (ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)756 zwz(ji,jj,jk) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 757 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 765 758 END DO 766 759 END DO … … 768 761 DO jj = 1, jpjm1 769 762 DO ji = 1, fs_jpim1 ! vector opt. 770 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv (ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) &771 & - e1u(ji ,jj+1) * pu (ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) &763 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 764 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 772 765 & * r1_e1e2f(ji,jj) ) 773 766 END DO … … 776 769 DO jj = 1, jpjm1 777 770 DO ji = 1, fs_jpim1 ! vector opt. 778 zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv (ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &779 & - ( pu (ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)771 zwz(ji,jj,jk) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 772 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 780 773 END DO 781 774 END DO … … 798 791 799 792 ! !== horizontal fluxes ==! 800 zwx(:,:) = e2u(:,:) * e3u (:,:,jk,ktlev) * pu(:,:,jk)801 zwy(:,:) = e1v(:,:) * e3v (:,:,jk,ktlev) * pv(:,:,jk)793 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 794 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 802 795 803 796 ! !== compute and add the vorticity term trend =! … … 805 798 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 806 799 DO ji = 2, jpi ! split in 2 parts due to vector opt. 807 z1_e3t = 1._wp / e3t (ji,jj,jk,ktlev)800 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 808 801 ztne(ji,jj) = ( zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) ) * z1_e3t 809 802 ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) ) * z1_e3t … … 813 806 DO jj = 3, jpj 814 807 DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3 815 z1_e3t = 1._wp / e3t (ji,jj,jk,ktlev)808 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 816 809 ztne(ji,jj) = ( zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) ) * z1_e3t 817 810 ztnw(ji,jj) = ( zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) ) * z1_e3t … … 826 819 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 827 820 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 828 pu _rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua829 pva _rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zva821 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 822 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 830 823 END DO 831 824 END DO … … 873 866 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 874 867 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 875 WRITE(numout,*) ' masked (=T) or u u(:,:,:,ktlev)masked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk868 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 876 869 ENDIF 877 870 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzad.F90
r10789 r10874 36 36 CONTAINS 37 37 38 SUBROUTINE dyn_zad ( kt , ktlev, pu_rhs, pv_rhs)38 SUBROUTINE dyn_zad ( kt ) 39 39 !!---------------------------------------------------------------------- 40 40 !! *** ROUTINE dynzad *** … … 44 44 !! 45 45 !! ** Method : The now vertical advection of momentum is given by: 46 !! w dz(u) = pu_rhs + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*ww) dk(uu(:,:,:,ktlev)) ]47 !! w dz(v) = pv_rhs + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*ww) dk(vv(:,:,:,ktlev)) ]48 !! Add this trend to the general trend ( pu_rhs,pv_rhs):49 !! ( pu_rhs,pv_rhs) = (pu_rhs,pv_rhs) + w dz(u,v)46 !! w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ] 47 !! w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ] 48 !! Add this trend to the general trend (ua,va): 49 !! (ua,va) = (ua,va) + w dz(u,v) 50 50 !! 51 !! ** Action : - Update ( pu_rhs,pv_rhs) with the vert. momentum adv. trends51 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 52 52 !! - Send the trends to trddyn for diagnostics (l_trddyn=T) 53 53 !!---------------------------------------------------------------------- 54 INTEGER, INTENT(in) :: kt ! ocean time-step index 55 INTEGER, INTENT(in) :: ktlev ! time level index for source terms 56 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends 54 INTEGER, INTENT(in) :: kt ! ocean time-step inedx 57 55 ! 58 56 INTEGER :: ji, jj, jk ! dummy loop indices … … 70 68 ENDIF 71 69 72 IF( l_trddyn ) THEN ! Save pu_rhs and pv_rhstrends70 IF( l_trddyn ) THEN ! Save ua and va trends 73 71 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 74 ztrdu(:,:,:) = pu_rhs(:,:,:)75 ztrdv(:,:,:) = pv_rhs(:,:,:)72 ztrdu(:,:,:) = ua(:,:,:) 73 ztrdv(:,:,:) = va(:,:,:) 76 74 ENDIF 77 75 … … 79 77 DO jj = 2, jpj ! vertical fluxes 80 78 DO ji = fs_2, jpi ! vector opt. 81 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * w w(ji,jj,jk)79 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 82 80 END DO 83 81 END DO 84 82 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 85 83 DO ji = fs_2, fs_jpim1 ! vector opt. 86 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( u u(ji,jj,jk-1,ktlev) - uu(ji,jj,jk,ktlev) )87 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( v v(ji,jj,jk-1,ktlev) - vv(ji,jj,jk,ktlev) )84 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) 85 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) 88 86 END DO 89 87 END DO … … 103 101 DO jj = 2, jpjm1 104 102 DO ji = fs_2, fs_jpim1 ! vector opt. 105 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,ktlev)106 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,ktlev)103 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 104 va(ji,jj,jk) = va(ji,jj,jk) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 107 105 END DO 108 106 END DO … … 110 108 111 109 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic 112 ztrdu(:,:,:) = pu_rhs(:,:,:) - ztrdu(:,:,:)113 ztrdv(:,:,:) = pv_rhs(:,:,:) - ztrdv(:,:,:)110 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 111 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 114 112 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 115 113 DEALLOCATE( ztrdu, ztrdv ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
r10825 r10874 45 45 CONTAINS 46 46 47 SUBROUTINE dyn_zdf( kt , ktlev1, ktlev2, ktlev3, kt2lev, pu_rhs, pv_rhs)47 SUBROUTINE dyn_zdf( kt ) 48 48 !!---------------------------------------------------------------------- 49 49 !! *** ROUTINE dyn_zdf *** … … 54 54 !! 55 55 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 !! pu_rhs = uu(:,:,:,ktlev1) + 2*dt * pu_rhsvector form or linear free surf.57 !! pu_rhs = ( e3u_b*uu(:,:,:,ktlev1) + 2*dt * e3u_n*pu_rhs ) / e3u(:,:,:,ktlev3)otherwise56 !! ua = ub + 2*dt * ua vector form or linear free surf. 57 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 58 58 !! - update the after velocity with the implicit vertical mixing. 59 59 !! This requires to solver the following system: 60 !! pu_rhs = pu_rhs + 1/e3u(:,:,:,ktlev3)dk+1[ mi(avm) / e3uw_a dk[ua] ]60 !! ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 61 61 !! with the following surface/top/bottom boundary condition: 62 62 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 63 63 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 64 64 !! 65 !! ** Action : ( pu_rhs,pv_rhs) after velocity65 !! ** Action : (ua,va) after velocity 66 66 !!--------------------------------------------------------------------- 67 INTEGER, INTENT(in) :: kt ! ocean time-step index 68 INTEGER, INTENT(in) :: ktlev1, ktlev2, ktlev3 ! time level indices for 3-time-level source terms 69 INTEGER, INTENT(in) :: kt2lev ! time level index for 2-time-level source terms 70 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends -> momentum after fields 67 INTEGER, INTENT(in) :: kt ! ocean time-step index 71 68 ! 72 69 INTEGER :: ji, jj, jk ! dummy loop indices … … 99 96 ! 100 97 ! !* explicit top/bottom drag case 101 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, u u(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs ) ! add top/bottom friction trend to (pu_rhs,pv_rhs)98 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, ub, vb, ua, va ) ! add top/bottom friction trend to (ua,va) 102 99 ! 103 100 ! 104 101 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 105 102 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 106 ztrdu(:,:,:) = pu_rhs(:,:,:)107 ztrdv(:,:,:) = pv_rhs(:,:,:)108 ENDIF 109 ! 110 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in pu_rhs,pv_rhs)103 ztrdu(:,:,:) = ua(:,:,:) 104 ztrdv(:,:,:) = va(:,:,:) 105 ENDIF 106 ! 107 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in ua,va) 111 108 ! 112 109 ! ! time stepping except vertical diffusion 113 110 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 114 111 DO jk = 1, jpkm1 115 pu_rhs(:,:,jk) = ( uu(:,:,jk,ktlev1) + r2dt * pu_rhs(:,:,jk) ) * umask(:,:,jk)116 pv_rhs(:,:,jk) = ( vv(:,:,jk,ktlev1) + r2dt * pv_rhs(:,:,jk) ) * vmask(:,:,jk)112 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 113 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 117 114 END DO 118 115 ELSE ! applied on thickness weighted velocity 119 116 DO jk = 1, jpkm1 120 pu_rhs(:,:,jk) = ( e3u(:,:,jk,ktlev1) * uu(:,:,jk,ktlev1) &121 & + r2dt * e3u (:,:,jk,ktlev2) * pu_rhs(:,:,jk) ) / e3u(:,:,jk,ktlev3) * umask(:,:,jk)122 pv_rhs(:,:,jk) = ( e3v(:,:,jk,ktlev1) * vv(:,:,jk,ktlev1) &123 & + r2dt * e3v (:,:,jk,ktlev2) * pv_rhs(:,:,jk) ) / e3v(:,:,jk,ktlev3) * vmask(:,:,jk)117 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 118 & + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 119 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 120 & + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 124 121 END DO 125 122 ENDIF … … 128 125 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 129 126 ! not lead to the effective stress seen over the whole barotropic loop. 130 ! G. Madec : in linear free surface, e3u (:,:,:,ktlev3) = e3u(:,:,:,ktlev2) = e3u_0, so systematic use of e3u(:,:,:,ktlev3)127 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 131 128 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 132 129 DO jk = 1, jpkm1 ! remove barotropic velocities 133 pu_rhs(:,:,jk) = ( pu_rhs(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)134 pv_rhs(:,:,jk) = ( pv_rhs(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)130 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 131 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 135 132 END DO 136 133 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only … … 138 135 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 139 136 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 140 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3)141 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3)142 pu_rhs(ji,jj,iku) = pu_rhs(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua143 pv_rhs(ji,jj,ikv) = pv_rhs(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va137 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 138 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 139 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 140 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 144 141 END DO 145 142 END DO … … 149 146 iku = miku(ji,jj) ! top ocean level at u- and v-points 150 147 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 151 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3)152 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3)153 pu_rhs(ji,jj,iku) = pu_rhs(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua154 pv_rhs(ji,jj,ikv) = pv_rhs(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va148 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 149 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 150 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 151 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 155 152 END DO 156 153 END DO … … 168 165 DO jj = 2, jpjm1 169 166 DO ji = fs_2, fs_jpim1 ! vector opt. 170 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point167 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 171 168 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 172 & / ( ze3ua * e3uw (ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk )169 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 173 170 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 174 & / ( ze3ua * e3uw (ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1)171 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 175 172 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 176 173 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 185 182 DO jj = 2, jpjm1 186 183 DO ji = fs_2, fs_jpim1 ! vector opt. 187 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point188 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw (ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk )189 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw (ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1)184 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 185 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 186 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 190 187 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 191 188 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 200 197 DO ji = fs_2, fs_jpim1 ! vector opt. 201 198 zwi(ji,jj,1) = 0._wp 202 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,1,ktlev2) + r_vvl * e3u(ji,jj,1,ktlev3)203 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw (ji,jj,2,kt2lev) ) * wumask(ji,jj,2)199 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 200 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 204 201 zWus = 0.5_wp * ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) 205 202 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) … … 213 210 DO jj = 2, jpjm1 214 211 DO ji = fs_2, fs_jpim1 ! vector opt. 215 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point212 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 216 213 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 217 & / ( ze3ua * e3uw (ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk )214 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 218 215 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 219 & / ( ze3ua * e3uw (ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1)216 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 220 217 zwi(ji,jj,jk) = zzwi 221 218 zws(ji,jj,jk) = zzws … … 228 225 DO jj = 2, jpjm1 229 226 DO ji = fs_2, fs_jpim1 ! vector opt. 230 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point231 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw (ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk )232 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw (ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1)227 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 228 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 229 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 233 230 zwi(ji,jj,jk) = zzwi 234 231 zws(ji,jj,jk) = zzws … … 257 254 DO ji = 2, jpim1 258 255 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 259 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) ! after scale factor at T-point256 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 260 257 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 261 258 END DO … … 266 263 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 267 264 iku = miku(ji,jj) ! ocean top level at u- and v-points 268 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) ! after scale factor at T-point265 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 269 266 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 270 267 END DO … … 285 282 ! m is decomposed in the product of an upper and a lower triangular matrix 286 283 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 287 ! The solution (the after velocity) is in pu_rhs284 ! The solution (the after velocity) is in ua 288 285 !----------------------------------------------------------------------- 289 286 ! … … 298 295 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 299 296 DO ji = fs_2, fs_jpim1 ! vector opt. 300 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,1,ktlev2) + r_vvl * e3u(ji,jj,1,ktlev3)301 pu_rhs(ji,jj,1) = pu_rhs(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &297 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 298 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 302 299 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 303 300 END DO … … 306 303 DO jj = 2, jpjm1 307 304 DO ji = fs_2, fs_jpim1 308 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pu_rhs(ji,jj,jk-1)305 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 309 306 END DO 310 307 END DO … … 313 310 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 314 311 DO ji = fs_2, fs_jpim1 ! vector opt. 315 pu_rhs(ji,jj,jpkm1) = pu_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)312 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 316 313 END DO 317 314 END DO … … 319 316 DO jj = 2, jpjm1 320 317 DO ji = fs_2, fs_jpim1 321 pu_rhs(ji,jj,jk) = ( pu_rhs(ji,jj,jk) - zws(ji,jj,jk) * pu_rhs(ji,jj,jk+1) ) / zwd(ji,jj,jk)318 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 322 319 END DO 323 320 END DO … … 334 331 DO jj = 2, jpjm1 335 332 DO ji = fs_2, fs_jpim1 ! vector opt. 336 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point333 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 337 334 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 338 & / ( ze3va * e3vw (ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk )335 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 339 336 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 340 & / ( ze3va * e3vw (ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1)337 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 341 338 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 342 339 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 351 348 DO jj = 2, jpjm1 352 349 DO ji = fs_2, fs_jpim1 ! vector opt. 353 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point354 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw (ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk )355 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw (ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1)350 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 351 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 352 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 356 353 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 357 354 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 366 363 DO ji = fs_2, fs_jpim1 ! vector opt. 367 364 zwi(ji,jj,1) = 0._wp 368 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,1,ktlev2) + r_vvl * e3v(ji,jj,1,ktlev3)369 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw (ji,jj,2,kt2lev) ) * wvmask(ji,jj,2)365 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 366 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 370 367 zWvs = 0.5_wp * ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) 371 368 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) … … 379 376 DO jj = 2, jpjm1 380 377 DO ji = fs_2, fs_jpim1 ! vector opt. 381 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point378 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 382 379 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 383 & / ( ze3va * e3vw (ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk )380 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 384 381 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 385 & / ( ze3va * e3vw (ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1)382 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 386 383 zwi(ji,jj,jk) = zzwi 387 384 zws(ji,jj,jk) = zzws … … 394 391 DO jj = 2, jpjm1 395 392 DO ji = fs_2, fs_jpim1 ! vector opt. 396 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point397 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw (ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk )398 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw (ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1)393 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 394 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 395 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 399 396 zwi(ji,jj,jk) = zzwi 400 397 zws(ji,jj,jk) = zzws … … 422 419 DO ji = 2, jpim1 423 420 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 424 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) ! after scale factor at T-point421 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 425 422 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 426 423 END DO … … 430 427 DO ji = 2, jpim1 431 428 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 432 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) ! after scale factor at T-point429 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 433 430 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 434 431 END DO … … 449 446 ! m is decomposed in the product of an upper and lower triangular matrix 450 447 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 451 ! The solution (after velocity) is in 2d array pv_rhs448 ! The solution (after velocity) is in 2d array va 452 449 !----------------------------------------------------------------------- 453 450 ! … … 462 459 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 463 460 DO ji = fs_2, fs_jpim1 ! vector opt. 464 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,1,ktlev2) + r_vvl * e3v(ji,jj,1,ktlev3)465 pv_rhs(ji,jj,1) = pv_rhs(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &461 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 462 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 466 463 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 467 464 END DO … … 470 467 DO jj = 2, jpjm1 471 468 DO ji = fs_2, fs_jpim1 ! vector opt. 472 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pv_rhs(ji,jj,jk-1)469 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 473 470 END DO 474 471 END DO … … 477 474 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 478 475 DO ji = fs_2, fs_jpim1 ! vector opt. 479 pv_rhs(ji,jj,jpkm1) = pv_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)476 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 480 477 END DO 481 478 END DO … … 483 480 DO jj = 2, jpjm1 484 481 DO ji = fs_2, fs_jpim1 485 pv_rhs(ji,jj,jk) = ( pv_rhs(ji,jj,jk) - zws(ji,jj,jk) * pv_rhs(ji,jj,jk+1) ) / zwd(ji,jj,jk)482 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 486 483 END DO 487 484 END DO … … 489 486 ! 490 487 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 491 ztrdu(:,:,:) = ( pu_rhs(:,:,:) - uu(:,:,:,ktlev1) ) / r2dt - ztrdu(:,:,:)492 ztrdv(:,:,:) = ( pv_rhs(:,:,:) - vv(:,:,:,ktlev1) ) / r2dt - ztrdv(:,:,:)488 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 489 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 493 490 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 494 491 DEALLOCATE( ztrdu, ztrdv )
Note: See TracChangeset
for help on using the changeset viewer.
