Changeset 5323
- Timestamp:
- 2015-06-01T09:59:54+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r5321 r5323 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, jpk,zhke )84 CALL wrk_alloc( jpi,jpj, 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(:,: ,jpk) = 0._wp99 97 98 zhke(:,:) = 0._wp 99 100 100 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 101 101 ! … … 108 108 zv = vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 109 109 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) 110 zhke(ji,jj ,jk) = 0.25_wp * ( zv + zu )111 END DO 110 zhke(ji,jj) = 0.25_wp * ( zv + zu ) 111 END DO 112 112 END DO 113 113 END DO … … 115 115 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 116 116 DO jk = 1, jpkm1 117 DO jj = 2, jpjm1 117 DO jj = 2, jpjm1 118 118 DO ji = fs_2, jpim1 ! vector opt. 119 119 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & … … 126 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) ) & 127 127 & + ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) * ( vn(ji-1,jj ,jk) + vn(ji+1,jj ,jk) ) 128 zhke(ji,jj ,jk) = r1_48 * ( zv + zu )129 END DO 128 zhke(ji,jj) = r1_48 * ( zv + zu ) 129 END DO 130 130 END DO 131 131 END DO … … 137 137 DO jj = 2, jpjm1 138 138 DO ji = fs_2, fs_jpim1 ! vector opt. 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 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) / e1u(ji,jj) 140 va(ji,jj,jk) = va(ji,jj,jk) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) / e2v(ji,jj) 141 END DO 142 142 END DO 143 143 END DO … … 153 153 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 154 154 ! 155 CALL wrk_dealloc( jpi,jpj, jpk,zhke )155 CALL wrk_dealloc( jpi,jpj, zhke ) 156 156 ! 157 157 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg')
Note: See TracChangeset
for help on using the changeset viewer.