- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynkeg.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/DYN/dynkeg.F90
r10996 r13463 36 36 37 37 !! * Substitutions 38 # include " vectopt_loop_substitute.h90"38 # include "do_loop_substitute.h90" 39 39 !!---------------------------------------------------------------------- 40 40 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 44 44 CONTAINS 45 45 46 SUBROUTINE dyn_keg( kt, kscheme )46 SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE dyn_keg *** … … 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 * un^2 + ((u n(j+1)+un(j-1))/2)^2 )60 !! + mj-1( 2 * vn^2 + ((v n(i+1)+vn(i-1))/2)^2 ) ]59 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 ) 60 !! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ] 61 61 !! 62 62 !! Take its horizontal gradient and add it to the general momentum 63 !! trend (ua,va).64 !! u a = ua- 1/e1u di[ zhke ]65 !! v a = va- 1/e2v dj[ zhke ]63 !! trend. 64 !! u(rhs) = u(rhs) - 1/e1u di[ zhke ] 65 !! v(rhs) = v(rhs) - 1/e2v dj[ zhke ] 66 66 !! 67 !! ** Action : - Update the ( ua, va) with the hor. ke gradient trend67 !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend 68 68 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 69 69 !! … … 71 71 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 72 72 !!---------------------------------------------------------------------- 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 73 INTEGER , INTENT( in ) :: kt ! ocean time-step index 74 INTEGER , INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 75 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 76 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 75 77 ! 76 INTEGER :: ji, jj, jk, jb ! dummy loop indices 77 INTEGER :: ifu, ifv, igrd, ib_bdy ! local integers 78 INTEGER :: ji, jj, jk ! dummy loop indices 78 79 REAL(wp) :: zu, zv ! local scalars 79 80 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke 80 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 81 REAL(wp) :: zweightu, zweightv82 82 !!---------------------------------------------------------------------- 83 83 ! … … 92 92 IF( l_trddyn ) THEN ! Save the input trends 93 93 ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) ) 94 ztrdu(:,:,:) = ua(:,:,:)95 ztrdv(:,:,:) = va(:,:,:)94 ztrdu(:,:,:) = puu(:,:,:,Krhs) 95 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 96 96 ENDIF 97 97 … … 101 101 ! 102 102 CASE ( nkeg_C2 ) !-- Standard scheme --! 103 DO jk = 1, jpkm1 104 DO jj = 2, jpj 105 DO ji = fs_2, jpi ! vector opt. 106 zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 107 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) 108 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 109 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) 110 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 111 END DO 112 END DO 113 END DO 114 ! 115 IF (ln_bdy) THEN 116 ! Maria Luneva & Fred Wobus: July-2016 117 ! compensate for lack of turbulent kinetic energy on liquid bdy points 118 DO ib_bdy = 1, nb_bdy 119 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 120 igrd = 1 ! compensating null velocity on the bdy 121 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 122 ji = idx_bdy(ib_bdy)%nbi(jb,igrd) ! maximum extent : from 2 to jpi-1 123 jj = idx_bdy(ib_bdy)%nbj(jb,igrd) ! maximum extent : from 2 to jpj-1 124 DO jk = 1, jpkm1 125 zhke(ji,jj,jk) = 0._wp 126 zweightu = umask(ji-1,jj ,jk) + umask(ji,jj,jk) 127 zweightv = vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) 128 zu = un(ji-1,jj ,jk) * un(ji-1,jj ,jk) + un(ji ,jj ,jk) * un(ji ,jj ,jk) 129 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) 130 IF( zweightu > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zu / (2._wp * zweightu) 131 IF( zweightv > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zv / (2._wp * zweightv) 132 END DO 133 END DO 134 END IF 135 CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy ) ! send 2 and recv jpi, jpj used in the computation of the speed tendencies 136 END DO 137 END IF 138 ! 103 DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 104 zu = puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) & 105 & + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) 106 zv = pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) & 107 & + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm) 108 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 109 END_3D 139 110 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 140 DO jk = 1, jpkm1 141 DO jj = 2, jpjm1 142 DO ji = fs_2, jpim1 ! vector opt. 143 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 144 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) & 145 & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) & 146 & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) 147 ! 148 zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 149 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) & 150 & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) & 151 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) 152 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 153 END DO 154 END DO 155 END DO 156 IF (ln_bdy) THEN 157 ! Maria Luneva & Fred Wobus: July-2016 158 ! compensate for lack of turbulent kinetic energy on liquid bdy points 159 DO ib_bdy = 1, nb_bdy 160 IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 161 igrd = 1 ! compensation null velocity on land at the bdy 162 DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 163 ji = idx_bdy(ib_bdy)%nbi(jb,igrd) ! maximum extent : from 2 to jpi-1 164 jj = idx_bdy(ib_bdy)%nbj(jb,igrd) ! maximum extent : from 2 to jpj-1 165 DO jk = 1, jpkm1 166 zhke(ji,jj,jk) = 0._wp 167 zweightu = 8._wp * ( umask(ji-1,jj ,jk) + umask(ji ,jj ,jk) ) & 168 & + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji ,jj-1,jk) + umask(ji ,jj+1,jk) ) 169 zweightv = 8._wp * ( vmask(ji ,jj-1,jk) + vmask(ji ,jj-1,jk) ) & 170 & + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj ,jk) + vmask(ji+1,jj ,jk) ) 171 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 172 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) & 173 & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) & 174 & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) 175 zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 176 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) & 177 & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) & 178 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) 179 IF( zweightu > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zu / ( 2._wp * zweightu ) 180 IF( zweightv > 0._wp ) zhke(ji,jj,jk) = zhke(ji,jj,jk) + zv / ( 2._wp * zweightv ) 181 END DO 182 END DO 183 END IF 184 END DO 185 END IF 186 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 111 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 112 zu = 8._wp * ( puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) & 113 & + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) ) & 114 & + ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) & 115 & + ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) * ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) 116 ! 117 zv = 8._wp * ( pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) & 118 & + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm) ) & 119 & + ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) & 120 & + ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) * ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) 121 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 122 END_3D 123 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1.0_wp ) 187 124 ! 188 125 END SELECT 189 126 ! 190 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! 191 DO jj = 2, jpjm1 192 DO ji = fs_2, fs_jpim1 ! vector opt. 193 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 194 va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 195 END DO 196 END DO 197 END DO 127 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 128 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 129 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 130 END_3D 198 131 ! 199 132 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 200 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)201 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)202 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt )133 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) 134 ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:) 135 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt, Kmm ) 203 136 DEALLOCATE( ztrdu , ztrdv ) 204 137 ENDIF 205 138 ! 206 IF( ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' keg - Ua: ', mask1=umask, &207 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )139 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg - Ua: ', mask1=umask, & 140 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 208 141 ! 209 142 IF( ln_timing ) CALL timing_stop('dyn_keg')
Note: See TracChangeset
for help on using the changeset viewer.
