Changeset 5328
- Timestamp:
- 2015-06-01T15:07:18+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r5324 r5328 7 7 !! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg 8 8 !! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module 9 !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option 9 !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option 10 10 !!---------------------------------------------------------------------- 11 11 12 12 !!---------------------------------------------------------------------- 13 13 !! dyn_keg : update the momentum trend with the horizontal tke … … 29 29 30 30 PUBLIC dyn_keg ! routine called by step module 31 31 32 32 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) 33 33 INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983 34 34 ! 35 35 REAL(wp) :: r1_48 = 1._wp / 48._wp !: =1/(4*2*6) 36 36 37 37 !! * Substitutions 38 38 # include "vectopt_loop_substitute.h90" 39 39 !!---------------------------------------------------------------------- 40 40 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 41 !! $Id$ 41 !! $Id$ 42 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 43 !!---------------------------------------------------------------------- … … 49 49 !! 50 50 !! ** Purpose : Compute the now momentum trend due to the horizontal 51 !! gradient of the horizontal kinetic energy and add it to the 51 !! gradient of the horizontal kinetic energy and add it to the 52 52 !! general momentum trend. 53 53 !! 54 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 55 !! conserve kinetic energy. Compute the now horizontal kinetic energy 54 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 55 !! conserve kinetic energy. Compute the now horizontal kinetic energy 56 56 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 57 57 !! * kscheme = nkeg_HW : Hollingsworth correction following … … 59 59 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((un(j+1)+un(j-1))/2)^2 ) 60 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 63 !! trend (ua,va). … … 72 72 !!---------------------------------------------------------------------- 73 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 75 75 ! 76 76 INTEGER :: ji, jj, jk ! dummy loop indices 77 77 REAL(wp) :: zu, zv ! temporary scalars 78 REAL(wp), POINTER, DIMENSION(:,: ):: zhke79 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 78 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 80 80 !!---------------------------------------------------------------------- 81 81 ! 82 82 IF( nn_timing == 1 ) CALL timing_start('dyn_keg') 83 83 ! 84 CALL wrk_alloc( jpi,jpj, zhke )84 CALL wrk_alloc( jpi,jpj,jpk, zhke ) 85 85 ! 86 86 IF( kt == nit000 ) THEN … … 92 92 IF( l_trddyn ) THEN ! Save ua and va trends 93 93 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 94 ztrdu(:,:,:) = ua(:,:,:) 95 ztrdv(:,:,:) = va(:,:,:) 94 ztrdu(:,:,:) = ua(:,:,:) 95 ztrdv(:,:,:) = va(:,:,:) 96 96 ENDIF 97 98 zhke(:,:) = 0._wp 99 100 DO jk = 1, jpkm1 101 102 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 103 ! 104 CASE ( nkeg_C2 ) !-- Standard scheme --! 105 97 98 zhke(:,:,jpk) = 0._wp 99 100 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 ! 102 CASE ( nkeg_C2 ) !-- Standard scheme --! 103 DO jk = 1, jpkm1 106 104 DO jj = 2, jpj 107 105 DO ji = fs_2, jpi ! vector opt. … … 110 108 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 111 109 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) 112 zhke(ji,jj ) = 0.25_wp * ( zv + zu )113 END DO 110 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 111 END DO 114 112 END DO 113 END DO 115 114 ! 116 117 118 DO jj = 2, jpjm1 115 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 116 DO jk = 1, jpkm1 117 DO jj = 2, jpjm1 119 118 DO ji = fs_2, jpim1 ! vector opt. 120 119 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & … … 127 126 & + ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) & 128 127 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) 129 zhke(ji,jj ) = r1_48 * ( zv + zu )130 END DO 128 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 129 END DO 131 130 END DO 132 CALL lbc_lnk( zhke, 'T', 1. ) 133 ! 134 END SELECT 131 END DO 132 CALL lbc_lnk( zhke, 'T', 1. ) 135 133 ! 134 END SELECT 135 ! 136 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! 136 137 DO jj = 2, jpjm1 137 138 DO ji = fs_2, fs_jpim1 ! vector opt. 138 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) / e1u(ji,jj)139 va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1 ) - zhke(ji,jj) ) / e2v(ji,jj)140 END DO 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 140 va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 141 END DO 141 142 END DO 142 143 143 END DO 144 144 ! … … 153 153 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 154 154 ! 155 CALL wrk_dealloc( jpi,jpj, zhke )155 CALL wrk_dealloc( jpi,jpj,jpk, zhke ) 156 156 ! 157 157 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg')
Note: See TracChangeset
for help on using the changeset viewer.