Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r3294 r5965 4 4 !! Ocean dynamics: kinetic energy gradient trend 5 5 !!====================================================================== 6 !! History : 1.0 ! 87-09 (P. Andrich, m.-a. Foujols) Original code 7 !! 7.0 ! 97-05 (G. Madec) Split dynber into dynkeg and dynhpg 8 !! 9.0 ! 02-07 (G. Madec) F90: Free form and module 6 !! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code 7 !! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg 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 10 !!---------------------------------------------------------------------- 10 11 … … 14 15 USE oce ! ocean dynamics and tracers 15 16 USE dom_oce ! ocean space and time domain 16 USE trdmod ! ocean dynamics trends 17 USE trdmod_oce ! ocean variables trends 17 USE trd_oce ! trends: ocean variables 18 USE trddyn ! trend manager: dynamics 19 ! 18 20 USE in_out_manager ! I/O manager 21 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 19 22 USE lib_mpp ! MPP library 20 23 USE prtctl ! Print control … … 27 30 PUBLIC dyn_keg ! routine called by step module 28 31 32 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) 33 INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983 34 ! 35 REAL(wp) :: r1_48 = 1._wp / 48._wp !: =1/(4*2*6) 36 29 37 !! * Substitutions 30 38 # include "vectopt_loop_substitute.h90" 31 39 !!---------------------------------------------------------------------- 32 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)40 !! NEMO/OPA 3.6 , NEMO Consortium (2015) 33 41 !! $Id$ 34 42 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 36 44 CONTAINS 37 45 38 SUBROUTINE dyn_keg( kt )46 SUBROUTINE dyn_keg( kt, kscheme ) 39 47 !!---------------------------------------------------------------------- 40 48 !! *** ROUTINE dyn_keg *** … … 44 52 !! general momentum trend. 45 53 !! 46 !! ** Method : 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 47 56 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 57 !! * kscheme = nkeg_HW : Hollingsworth correction following 58 !! Arakawa (2001). The now horizontal kinetic energy is given by: 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 !! 48 62 !! Take its horizontal gradient and add it to the general momentum 49 63 !! trend (ua,va). … … 52 66 !! 53 67 !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 54 !! - save this trends (l_trddyn=T) for post-processing 68 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 69 !! 70 !! ** References : Arakawa, A., International Geophysics 2001. 71 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 55 72 !!---------------------------------------------------------------------- 56 INTEGER, INTENT( in ) :: kt ! ocean time-step index 57 !! 73 INTEGER, INTENT( in ) :: kt ! ocean time-step index 74 INTEGER, INTENT( in ) :: kscheme ! =0/1 type of KEG scheme 75 ! 58 76 INTEGER :: ji, jj, jk ! dummy loop indices 59 77 REAL(wp) :: zu, zv ! temporary scalars … … 62 80 !!---------------------------------------------------------------------- 63 81 ! 64 IF( nn_timing == 1 ) CALL timing_start('dyn_keg')82 IF( nn_timing == 1 ) CALL timing_start('dyn_keg') 65 83 ! 66 CALL wrk_alloc( jpi, jpj, jpk,zhke )84 CALL wrk_alloc( jpi,jpj,jpk, zhke ) 67 85 ! 68 86 IF( kt == nit000 ) THEN 69 87 IF(lwp) WRITE(numout,*) 70 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend '88 IF(lwp) WRITE(numout,*) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme 71 89 IF(lwp) WRITE(numout,*) '~~~~~~~' 72 90 ENDIF 73 91 74 92 IF( l_trddyn ) THEN ! Save ua and va trends 75 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )93 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 76 94 ztrdu(:,:,:) = ua(:,:,:) 77 95 ztrdv(:,:,:) = va(:,:,:) 78 96 ENDIF 79 97 80 ! ! =============== 81 DO jk = 1, jpkm1 ! Horizontal slab 82 ! ! =============== 83 DO jj = 2, jpj ! Horizontal kinetic energy at T-point 84 DO ji = fs_2, jpi ! vector opt. 85 zu = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 86 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) 87 zv = 0.25 * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 88 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 89 zhke(ji,jj,jk) = zv + zu 90 !!gm simplier coding ==>> ~ faster 91 ! don't forget to suppress local zu zv scalars 92 ! zhke(ji,jj,jk) = 0.25 * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 93 ! & + un(ji ,jj ,jk) * un(ji ,jj ,jk) & 94 ! & + vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 95 ! & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) 96 !!gm end <<== 97 END DO 98 END DO 99 DO jj = 2, jpjm1 ! add the gradient of kinetic energy to the general momentum trends 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 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 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 116 DO jk = 1, jpkm1 117 DO jj = 2, jpjm1 118 DO ji = fs_2, jpim1 ! vector opt. 119 zu = 8._wp * ( un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 120 & + un(ji ,jj ,jk) * un(ji ,jj ,jk) ) & 121 & + ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) & 122 & + ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) * ( un(ji ,jj-1,jk) + un(ji ,jj+1,jk) ) 123 ! 124 zv = 8._wp * ( vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 125 & + vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) & 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 & + ( 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 130 END DO 131 END DO 132 CALL lbc_lnk( zhke, 'T', 1. ) 133 ! 134 END SELECT 135 ! 136 DO jk = 1, jpkm1 !== grad( KE ) added to the general momentum trends ==! 137 DO jj = 2, jpjm1 100 138 DO ji = fs_2, fs_jpim1 ! vector opt. 101 139 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) … … 103 141 END DO 104 142 END DO 105 !!gm idea to be tested ==>> is it faster on scalar computers ? 106 ! DO jj = 2, jpjm1 ! add the gradient of kinetic energy to the general momentum trends 107 ! DO ji = fs_2, fs_jpim1 ! vector opt. 108 ! ua(ji,jj,jk) = ua(ji,jj,jk) - 0.25 * ( + un(ji+1,jj ,jk) * un(ji+1,jj ,jk) & 109 ! & + vn(ji+1,jj-1,jk) * vn(ji+1,jj-1,jk) & 110 ! & + vn(ji+1,jj ,jk) * vn(ji+1,jj ,jk) & 111 ! ! 112 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 113 ! & - vn(ji ,jj-1,jk) * vn(ji ,jj-1,jk) & 114 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e1u(ji,jj) 115 ! ! 116 ! va(ji,jj,jk) = va(ji,jj,jk) - 0.25 * ( un(ji-1,jj+1,jk) * un(ji-1,jj+1,jk) & 117 ! & + un(ji ,jj+1,jk) * un(ji ,jj+1,jk) & 118 ! & + vn(ji ,jj+1,jk) * vn(ji ,jj+1,jk) & 119 ! ! 120 ! & - un(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 121 ! & - un(ji ,jj ,jk) * un(ji ,jj ,jk) & 122 ! & - vn(ji ,jj ,jk) * vn(ji ,jj ,jk) ) / e2v(ji,jj) 123 ! END DO 124 ! END DO 125 !!gm en idea <<== 126 ! ! =============== 127 END DO ! End of slab 128 ! ! =============== 129 130 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 143 END DO 144 ! 145 IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic 131 146 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 132 147 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 133 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )134 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )148 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 149 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 135 150 ENDIF 136 151 ! … … 138 153 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 139 154 ! 140 CALL wrk_dealloc( jpi, jpj, jpk,zhke )155 CALL wrk_dealloc( jpi,jpj,jpk, zhke ) 141 156 ! 142 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg')157 IF( nn_timing == 1 ) CALL timing_stop('dyn_keg') 143 158 ! 144 159 END SUBROUTINE dyn_keg
Note: See TracChangeset
for help on using the changeset viewer.
