Changeset 5836 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2015-10-26T15:49:40+01:00 (8 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 5 deleted
- 16 edited
- 2 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r5322 r5836 76 76 CASE ( 3 ) 77 77 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 78 !79 CASE (-1 ) ! esopa: test all possibility with control print80 CALL dyn_keg ( kt, nn_dynkeg )81 CALL dyn_zad ( kt )82 CALL dyn_adv_cen2( kt )83 CALL dyn_adv_ubs ( kt )84 78 END SELECT 85 79 ! … … 126 120 IF( ln_dynadv_cen2 ) ioptio = ioptio + 1 127 121 IF( ln_dynadv_ubs ) ioptio = ioptio + 1 128 IF( lk_esopa ) ioptio = 1129 122 130 123 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) … … 139 132 IF( ln_dynadv_cen2 ) nadv = 2 140 133 IF( ln_dynadv_ubs ) nadv = 3 141 IF( lk_esopa ) nadv = -1142 134 143 135 IF(lwp) THEN ! Print the choice … … 151 143 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 152 144 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used' 153 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection formulation'154 145 ENDIF 155 146 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
r4990 r5836 90 90 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 91 91 DO ji = fs_2, fs_jpim1 ! vector opt. 92 zbu = e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk)93 zbv = e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk)92 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 93 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 94 94 ! 95 95 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & … … 114 114 DO jk = 1, jpkm1 ! ==================== ! 115 115 ! ! Vertical volume fluxesÊ 116 zfw(:,:,jk) = 0.25 * e1 t(:,:) *e2t(:,:) * wn(:,:,jk)116 zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk) 117 117 ! 118 118 IF( jk == 1 ) THEN ! surface/bottom advective fluxes … … 144 144 DO ji = fs_2, fs_jpim1 ! vector opt. 145 145 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 146 & / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )146 & / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 147 147 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 148 & / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )148 & / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 149 149 END DO 150 150 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r5069 r5836 181 181 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 182 182 DO ji = fs_2, fs_jpim1 ! vector opt. 183 zbu = e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk)184 zbv = e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk)183 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 184 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 185 185 ! 186 186 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & … … 203 203 DO jk = 1, jpkm1 ! ==================== ! 204 204 ! ! Vertical volume fluxesÊ 205 zfw(:,:,jk) = 0.25 * e1 t(:,:) *e2t(:,:) * wn(:,:,jk)205 zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk) 206 206 ! 207 207 IF( jk == 1 ) THEN ! surface/bottom advective fluxes … … 233 233 DO ji = fs_2, fs_jpim1 ! vector opt. 234 234 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 235 & / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )235 & / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 236 236 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 237 & / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )237 & / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 238 238 END DO 239 239 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5224 r5836 1046 1046 DO jj = 2, jpjm1 1047 1047 DO ji = 2, jpim1 1048 zsshu_n(ji,jj) = (e1 2u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * &1049 & r1_e1 2u(ji,jj) * umask(ji,jj,1) * 0.5_wp1050 zsshv_n(ji,jj) = (e1 2v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * &1051 & r1_e1 2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp1048 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1049 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1050 zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * & 1051 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1052 1052 END DO 1053 1053 END DO -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r4990 r5836 4 4 !! Ocean physics: lateral diffusivity trends 5 5 !!===================================================================== 6 !! History : 9.0 ! 05-11 (G. Madec) Original code (new step architecture) 6 !! History : 2.0 ! 2005-11 (G. Madec) Original code (new step architecture) 7 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 8 !! ! add velocity dependent coefficient and optional read in file 7 9 !!---------------------------------------------------------------------- 8 10 … … 14 16 USE dom_oce ! ocean space and time domain 15 17 USE phycst ! physical constants 16 USE ldfdyn_oce ! ocean dynamics lateral physics 17 USE ldftra_oce ! ocean tracers lateral physics 18 USE ldfslp ! lateral mixing: slopes of mixing orientation 19 USE dynldf_bilapg ! lateral mixing (dyn_ldf_bilapg routine) 20 USE dynldf_bilap ! lateral mixing (dyn_ldf_bilap routine) 21 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine) 22 USE dynldf_lap ! lateral mixing (dyn_ldf_lap routine) 18 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 19 USE ldfslp ! lateral diffusion: slopes of mixing orientation 20 USE dynldf_lap_blp ! lateral mixing (dyn_ldf_lap & dyn_ldf_blp routines) 21 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine ) 23 22 USE trd_oce ! trends: ocean variables 24 USE trddyn ! trend manager: dynamics (trd_dyn 23 USE trddyn ! trend manager: dynamics (trd_dyn routine) 25 24 ! 26 25 USE prtctl ! Print control … … 28 27 USE lib_mpp ! distribued memory computing library 29 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 30 USE wrk_nemo 31 USE timing 29 USE wrk_nemo ! Memory Allocation 30 USE timing ! Timing 32 31 33 32 IMPLICIT NONE … … 37 36 PUBLIC dyn_ldf_init ! called by opa module 38 37 39 INTEGER :: nldf = -2 ! type of lateral diffusion used defined from ln_dynldf_... namlist logicals) 38 ! ! Flag to control the type of lateral viscous operator 39 INTEGER, PARAMETER, PUBLIC :: np_ERROR =-10 ! error in setting the operator 40 INTEGER, PARAMETER, PUBLIC :: np_no_ldf = 00 ! without operator (i.e. no lateral viscous trend) 41 ! !! laplacian ! bilaplacian ! 42 INTEGER, PARAMETER, PUBLIC :: np_lap = 10 , np_blp = 20 ! iso-level operator 43 INTEGER, PARAMETER, PUBLIC :: np_lap_i = 11 ! iso-neutral or geopotential operator 44 45 INTEGER :: nldf ! type of lateral diffusion used defined from ln_dynldf_... (namlist logicals) 40 46 41 47 !! * Substitutions … … 43 49 # include "vectopt_loop_substitute.h90" 44 50 !!---------------------------------------------------------------------- 45 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)51 !! NEMO/OPA 3.7 , NEMO Consortium (2015) 46 52 !! $Id$ 47 53 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 62 68 IF( nn_timing == 1 ) CALL timing_start('dyn_ldf') 63 69 ! 64 IF( l_trddyn ) THEN ! temporary save of ta and satrends65 CALL wrk_alloc( jpi, jpj, jpk,ztrdu, ztrdv )70 IF( l_trddyn ) THEN ! temporary save of momentum trends 71 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 66 72 ztrdu(:,:,:) = ua(:,:,:) 67 73 ztrdv(:,:,:) = va(:,:,:) … … 70 76 SELECT CASE ( nldf ) ! compute lateral mixing trend and add it to the general trend 71 77 ! 72 CASE ( 0 ) ; CALL dyn_ldf_lap ( kt ) ! iso-level laplacian 73 CASE ( 1 ) ; CALL dyn_ldf_iso ( kt ) ! rotated laplacian (except dk[ dk[.] ] part) 74 CASE ( 2 ) ; CALL dyn_ldf_bilap ( kt ) ! iso-level bilaplacian 75 CASE ( 3 ) ; CALL dyn_ldf_bilapg ( kt ) ! s-coord. horizontal bilaplacian 76 CASE ( 4 ) ! iso-level laplacian + bilaplacian 77 CALL dyn_ldf_lap ( kt ) 78 CALL dyn_ldf_bilap ( kt ) 79 CASE ( 5 ) ! rotated laplacian + bilaplacian (s-coord) 80 CALL dyn_ldf_iso ( kt ) 81 CALL dyn_ldf_bilapg ( kt ) 78 CASE ( np_lap ) ; CALL dyn_ldf_lap ( kt, ub, vb, ua, va, 1 ) ! iso-level laplacian 79 CASE ( np_lap_i ) ; CALL dyn_ldf_iso ( kt ) ! rotated laplacian 80 CASE ( np_blp ) ; CALL dyn_ldf_blp ( kt, ub, vb, ua, va ) ! iso-level bi-laplacian 82 81 ! 83 CASE ( -1 ) ! esopa: test all possibility with control print84 CALL dyn_ldf_lap ( kt )85 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf0 - Ua: ', mask1=umask, &86 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )87 CALL dyn_ldf_iso ( kt )88 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf1 - Ua: ', mask1=umask, &89 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )90 CALL dyn_ldf_bilap ( kt )91 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf2 - Ua: ', mask1=umask, &92 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )93 CALL dyn_ldf_bilapg ( kt )94 CALL prt_ctl( tab3d_1=ua, clinfo1=' ldf3 - Ua: ', mask1=umask, &95 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )96 !97 CASE ( -2 ) ! neither laplacian nor bilaplacian schemes used98 IF( kt == nit000 ) THEN99 IF(lwp) WRITE(numout,*)100 IF(lwp) WRITE(numout,*) 'dyn_ldf : no lateral diffusion on momentum setup'101 IF(lwp) WRITE(numout,*) '~~~~~~~ '102 ENDIF103 82 END SELECT 104 83 … … 107 86 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 108 87 CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 109 CALL wrk_dealloc( jpi, jpj, jpk,ztrdu, ztrdv )88 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 110 89 ENDIF 111 90 ! ! print sum trends (used for debugging) … … 126 105 INTEGER :: ioptio, ierr ! temporary integers 127 106 !!---------------------------------------------------------------------- 128 107 ! 129 108 ! ! Namelist nam_dynldf: already read in ldfdyn module 130 109 ! 131 110 IF(lwp) THEN ! Namelist print 132 111 WRITE(numout,*) … … 134 113 WRITE(numout,*) '~~~~~~~~~~~' 135 114 WRITE(numout,*) ' Namelist nam_dynldf : set lateral mixing parameters (type, direction, coefficients)' 136 WRITE(numout,*) ' laplacian operator ln_dynldf_lap 137 WRITE(numout,*) ' bilaplacian operator ln_dynldf_b ilap = ', ln_dynldf_bilap138 WRITE(numout,*) ' iso-level ln_dynldf_lev el = ', ln_dynldf_level139 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor 140 WRITE(numout,*) ' iso-neutral ln_dynldf_iso 115 WRITE(numout,*) ' laplacian operator ln_dynldf_lap = ', ln_dynldf_lap 116 WRITE(numout,*) ' bilaplacian operator ln_dynldf_blp = ', ln_dynldf_blp 117 WRITE(numout,*) ' iso-level ln_dynldf_lev = ', ln_dynldf_lev 118 WRITE(numout,*) ' horizontal (geopotential) ln_dynldf_hor = ', ln_dynldf_hor 119 WRITE(numout,*) ' iso-neutral ln_dynldf_iso = ', ln_dynldf_iso 141 120 ENDIF 142 143 ! ! control the consistency121 ! ! use of lateral operator or not 122 nldf = np_ERROR 144 123 ioptio = 0 145 IF( ln_dynldf_lap ) ioptio = ioptio + 1 146 IF( ln_dynldf_bilap ) ioptio = ioptio + 1 147 IF( ioptio < 1 ) CALL ctl_warn( ' neither laplacian nor bilaplacian operator set for dynamics' ) 148 ioptio = 0 149 IF( ln_dynldf_level ) ioptio = ioptio + 1 150 IF( ln_dynldf_hor ) ioptio = ioptio + 1 151 IF( ln_dynldf_iso ) ioptio = ioptio + 1 152 IF( ioptio > 1 ) CALL ctl_stop( ' use only ONE direction (level/hor/iso)' ) 153 154 IF( ln_dynldf_iso .AND. ln_traldf_hor ) CALL ctl_stop & 155 & ( 'Not sensible to use geopotential diffusion for tracers with isoneutral diffusion for dynamics' ) 156 157 ! ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 158 ierr = 0 159 IF ( ln_dynldf_lap ) THEN ! laplacian operator 160 IF ( ln_zco ) THEN ! z-coordinate 161 IF ( ln_dynldf_level ) nldf = 0 ! iso-level (no rotation) 162 IF ( ln_dynldf_hor ) nldf = 0 ! horizontal (no rotation) 163 IF ( ln_dynldf_iso ) nldf = 1 ! isoneutral ( rotation) 124 IF( ln_dynldf_lap ) ioptio = ioptio + 1 125 IF( ln_dynldf_blp ) ioptio = ioptio + 1 126 IF( ioptio > 1 ) CALL ctl_stop( 'dyn_ldf_init: use ONE or NONE of the 2 lap/bilap operator type on momentum' ) 127 IF( ioptio == 0 ) nldf = np_no_ldf ! No lateral mixing operator 128 ! 129 IF( nldf /= np_no_ldf ) THEN ! direction ==>> type of operator 130 ioptio = 0 131 IF( ln_dynldf_lev ) ioptio = ioptio + 1 132 IF( ln_dynldf_hor ) ioptio = ioptio + 1 133 IF( ln_dynldf_iso ) ioptio = ioptio + 1 134 IF( ioptio > 1 ) CALL ctl_stop( ' use only ONE direction (level/hor/iso)' ) 135 IF( ioptio == 0 ) CALL ctl_stop( ' use at least ONE direction (level/hor/iso)' ) 136 ! 137 ! ! Set nldf, the type of lateral diffusion, from ln_dynldf_... logicals 138 ierr = 0 139 IF ( ln_dynldf_lap ) THEN ! laplacian operator 140 IF ( ln_zco ) THEN ! z-coordinate 141 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 142 IF ( ln_dynldf_hor ) nldf = np_lap ! iso-level = horizontal (no rotation) 143 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 144 ENDIF 145 IF ( ln_zps ) THEN ! z-coordinate with partial step 146 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level (no rotation) 147 IF ( ln_dynldf_hor ) nldf = np_lap ! iso-level (no rotation) 148 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 149 ENDIF 150 IF ( ln_sco ) THEN ! s-coordinate 151 IF ( ln_dynldf_lev ) nldf = np_lap ! iso-level = horizontal (no rotation) 152 IF ( ln_dynldf_hor ) nldf = np_lap_i ! horizontal ( rotation) 153 IF ( ln_dynldf_iso ) nldf = np_lap_i ! iso-neutral ( rotation) 154 ENDIF 164 155 ENDIF 165 IF ( ln_zps ) THEN ! z-coordinate 166 IF ( ln_dynldf_level ) ierr = 1 ! iso-level not allowed 167 IF ( ln_dynldf_hor ) nldf = 0 ! horizontal (no rotation) 168 IF ( ln_dynldf_iso ) nldf = 1 ! isoneutral ( rotation) 156 ! 157 IF( ln_dynldf_blp ) THEN ! bilaplacian operator 158 IF ( ln_zco ) THEN ! z-coordinate 159 IF ( ln_dynldf_lev ) nldf = np_blp ! iso-level = horizontal (no rotation) 160 IF ( ln_dynldf_hor ) nldf = np_blp ! iso-level = horizontal (no rotation) 161 IF ( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 162 ENDIF 163 IF ( ln_zps ) THEN ! z-coordinate with partial step 164 IF ( ln_dynldf_lev ) nldf = np_blp ! iso-level (no rotation) 165 IF ( ln_dynldf_hor ) nldf = np_blp ! iso-level (no rotation) 166 IF ( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 167 ENDIF 168 IF ( ln_sco ) THEN ! s-coordinate 169 IF ( ln_dynldf_lev ) nldf = np_blp ! iso-level (no rotation) 170 IF ( ln_dynldf_hor ) ierr = 2 ! horizontal ( rotation) 171 IF ( ln_dynldf_iso ) ierr = 2 ! iso-neutral ( rotation) 172 ENDIF 169 173 ENDIF 170 IF ( ln_sco ) THEN ! s-coordinate 171 IF ( ln_dynldf_level ) nldf = 0 ! iso-level (no rotation) 172 IF ( ln_dynldf_hor ) nldf = 1 ! horizontal ( rotation) 173 IF ( ln_dynldf_iso ) nldf = 1 ! isoneutral ( rotation) 174 ENDIF 175 ENDIF 176 177 IF( ln_dynldf_bilap ) THEN ! bilaplacian operator 178 IF ( ln_zco ) THEN ! z-coordinate 179 IF ( ln_dynldf_level ) nldf = 2 ! iso-level (no rotation) 180 IF ( ln_dynldf_hor ) nldf = 2 ! horizontal (no rotation) 181 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 182 ENDIF 183 IF ( ln_zps ) THEN ! z-coordinate 184 IF ( ln_dynldf_level ) ierr = 1 ! iso-level not allowed 185 IF ( ln_dynldf_hor ) nldf = 2 ! horizontal (no rotation) 186 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 187 ENDIF 188 IF ( ln_sco ) THEN ! s-coordinate 189 IF ( ln_dynldf_level ) nldf = 2 ! iso-level (no rotation) 190 IF ( ln_dynldf_hor ) nldf = 3 ! horizontal ( rotation) 191 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 192 ENDIF 193 ENDIF 194 195 IF( ln_dynldf_lap .AND. ln_dynldf_bilap ) THEN ! mixed laplacian and bilaplacian operators 196 IF ( ln_zco ) THEN ! z-coordinate 197 IF ( ln_dynldf_level ) nldf = 4 ! iso-level (no rotation) 198 IF ( ln_dynldf_hor ) nldf = 4 ! horizontal (no rotation) 199 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 200 ENDIF 201 IF ( ln_zps ) THEN ! z-coordinate 202 IF ( ln_dynldf_level ) ierr = 1 ! iso-level not allowed 203 IF ( ln_dynldf_hor ) nldf = 4 ! horizontal (no rotation) 204 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 205 ENDIF 206 IF ( ln_sco ) THEN ! s-coordinate 207 IF ( ln_dynldf_level ) nldf = 4 ! iso-level (no rotation) 208 IF ( ln_dynldf_hor ) nldf = 5 ! horizontal ( rotation) 209 IF ( ln_dynldf_iso ) ierr = 2 ! isoneutral ( rotation) 210 ENDIF 211 ENDIF 212 213 IF( lk_esopa ) nldf = -1 ! esopa test 214 215 IF( ierr == 1 ) CALL ctl_stop( 'iso-level in z-coordinate - partial step, not allowed' ) 216 IF( ierr == 2 ) CALL ctl_stop( 'isoneutral bilaplacian operator does not exist' ) 217 IF( nldf == 1 .OR. nldf == 3 ) THEN ! rotation 218 IF( .NOT.lk_ldfslp ) CALL ctl_stop( 'the rotation of the diffusive tensor require key_ldfslp' ) 174 ! 175 IF( ierr == 2 ) CALL ctl_stop( 'rotated bi-laplacian operator does not exist' ) 176 ! 177 IF( nldf == np_lap_i ) l_ldfslp = .TRUE. ! rotation require the computation of the slopes 178 ! 219 179 ENDIF 220 180 221 181 IF(lwp) THEN 222 182 WRITE(numout,*) 223 IF( nldf == -2 ) WRITE(numout,*) ' neither laplacian nor bilaplacian schemes used' 224 IF( nldf == -1 ) WRITE(numout,*) ' ESOPA test All scheme used' 225 IF( nldf == 0 ) WRITE(numout,*) ' laplacian operator' 226 IF( nldf == 1 ) WRITE(numout,*) ' rotated laplacian operator' 227 IF( nldf == 2 ) WRITE(numout,*) ' bilaplacian operator' 228 IF( nldf == 3 ) WRITE(numout,*) ' rotated bilaplacian' 229 IF( nldf == 4 ) WRITE(numout,*) ' laplacian and bilaplacian operators' 230 IF( nldf == 5 ) WRITE(numout,*) ' rotated laplacian and bilaplacian operators' 183 IF( nldf == np_no_ldf ) WRITE(numout,*) ' NO lateral viscosity' 184 IF( nldf == np_lap ) WRITE(numout,*) ' iso-level laplacian operator' 185 IF( nldf == np_lap_i ) WRITE(numout,*) ' rotated laplacian operator with iso-level background' 186 IF( nldf == np_blp ) WRITE(numout,*) ' iso-level bi-laplacian operator' 231 187 ENDIF 232 188 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r4990 r5836 2 2 !!====================================================================== 3 3 !! *** MODULE dynldf_iso *** 4 !! Ocean dynamics: lateral viscosity trend4 !! Ocean dynamics: lateral viscosity trend (rotated laplacian operator) 5 5 !!====================================================================== 6 6 !! History : OPA ! 97-07 (G. Madec) Original code … … 8 8 !! - ! 2004-08 (C. Talandier) New trends organization 9 9 !! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion 10 !! 3.7 ! 2014-01 (F. Lemarie, G. Madec) restructuration/simplification of ahm specification, 11 !! ! add velocity dependent coefficient and optional read in file 10 12 !!---------------------------------------------------------------------- 11 #if defined key_ldfslp || defined key_esopa 12 !!---------------------------------------------------------------------- 13 !! 'key_ldfslp' slopes of the direction of mixing 13 14 14 !!---------------------------------------------------------------------- 15 15 !! dyn_ldf_iso : update the momentum trend with the horizontal part … … 19 19 USE oce ! ocean dynamics and tracers 20 20 USE dom_oce ! ocean space and time domain 21 USE ldfdyn _oce ! ocean dynamics lateral physics22 USE ldftra _oce ! ocean tracer lateral physics21 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 22 USE ldftra ! lateral physics: eddy diffusivity 23 23 USE zdf_oce ! ocean vertical physics 24 24 USE ldfslp ! iso-neutral slopes 25 25 ! 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link)27 26 USE in_out_manager ! I/O manager 28 27 USE lib_mpp ! MPP library 28 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 29 29 USE prtctl ! Print control 30 30 USE wrk_nemo ! Memory Allocation … … 42 42 !! * Substitutions 43 43 # include "domzgr_substitute.h90" 44 # include "ldfdyn_substitute.h90"45 44 # include "vectopt_loop_substitute.h90" 46 45 !!---------------------------------------------------------------------- … … 83 82 !! horizontal fluxes associated with the rotated lateral mixing: 84 83 !! u-component: 85 !! ziut = ( ahmt + ahmb0) e2t * e3t / e1t di[ ub ]86 !! - ahmte2t * mi-1(uslp) dk[ mi(mk(ub)) ]87 !! zjuf = ( ahmf + ahmb0) e1f * e3f / e2f dj[ ub ]88 !! - ahmfe1f * mi(vslp) dk[ mj(mk(ub)) ]84 !! ziut = ( ahmt + rn_ahm_b ) e2t * e3t / e1t di[ ub ] 85 !! - ahmt e2t * mi-1(uslp) dk[ mi(mk(ub)) ] 86 !! zjuf = ( ahmf + rn_ahm_b ) e1f * e3f / e2f dj[ ub ] 87 !! - ahmf e1f * mi(vslp) dk[ mj(mk(ub)) ] 89 88 !! v-component: 90 !! zivf = ( ahmf + ahmb0) e2t * e3t / e1t di[ vb ]91 !! - ahmfe2t * mj(uslp) dk[ mi(mk(vb)) ]92 !! zjvt = ( ahmt + ahmb0) e1f * e3f / e2f dj[ ub ]93 !! - ahmte1f * mj-1(vslp) dk[ mj(mk(vb)) ]89 !! zivf = ( ahmf + rn_ahm_b ) e2t * e3t / e1t di[ vb ] 90 !! - ahmf e2t * mj(uslp) dk[ mi(mk(vb)) ] 91 !! zjvt = ( ahmt + rn_ahm_b ) e1f * e3f / e2f dj[ ub ] 92 !! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vb)) ] 94 93 !! take the horizontal divergence of the fluxes: 95 94 !! diffu = 1/(e1u*e2u*e3u) { di [ ziut ] + dj-1[ zjuf ] } … … 106 105 !! of the rotated operator in dynzdf module 107 106 !!---------------------------------------------------------------------- 108 !109 107 INTEGER, INTENT( in ) :: kt ! ocean time-step index 110 108 ! 111 109 INTEGER :: ji, jj, jk ! dummy loop indices 112 110 REAL(wp) :: zabe1, zabe2, zcof1, zcof2 ! local scalars 113 REAL(wp) :: zmskt, zmskf , zbu, zbv, zuah, zvah! - -111 REAL(wp) :: zmskt, zmskf ! - - 114 112 REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - - 115 113 REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - - … … 130 128 ENDIF 131 129 132 ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum 130 !!gm bug is dyn_ldf_iso called before tra_ldf_iso .... <<<<<===== TO BE CHECKED 131 ! s-coordinate: Iso-level diffusion on momentum but not on tracer 133 132 IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN 134 133 ! … … 136 135 DO jj = 2, jpjm1 137 136 DO ji = 2, jpim1 138 uslp (ji,jj,jk) = - 1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk)) * umask(ji,jj,jk)139 vslp (ji,jj,jk) = - 1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk)) * vmask(ji,jj,jk)140 wslpi(ji,jj,jk) = - 1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk)) * tmask(ji,jj,jk) * 0.5141 wslpj(ji,jj,jk) = - 1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk)) * tmask(ji,jj,jk) * 0.5137 uslp (ji,jj,jk) = - ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 138 vslp (ji,jj,jk) = - ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 139 wslpi(ji,jj,jk) = - ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * r1_e1t(ji,jj) * tmask(ji,jj,jk) * 0.5 140 wslpj(ji,jj,jk) = - ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * r1_e2t(ji,jj) * tmask(ji,jj,jk) * 0.5 142 141 END DO 143 142 END DO … … 184 183 DO jj = 2, jpjm1 185 184 DO ji = fs_2, jpi ! vector opt. 186 zabe1 = ( fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) /e1t(ji,jj)187 188 zmskt = 1. /MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1)&189 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1.)190 191 zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )185 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) * r1_e1t(ji,jj) 186 187 zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & 188 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ) , 1._wp ) 189 190 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 192 191 192 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 193 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & 194 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk) 195 END DO 196 END DO 197 ELSE ! other coordinate system (zco or sco) : e3t 198 DO jj = 2, jpjm1 199 DO ji = fs_2, jpi ! vector opt. 200 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * fse3t(ji,jj,jk) * r1_e1t(ji,jj) 201 202 zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) & 203 & + umask(ji-1,jj,jk+1) + umask(ji,jj,jk ) , 1._wp ) 204 205 zcof1 = - rn_aht_0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) 206 193 207 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) & 194 208 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) & … … 196 210 END DO 197 211 END DO 198 ELSE ! other coordinate system (zco or sco) : e3t199 DO jj = 2, jpjm1200 DO ji = fs_2, jpi ! vector opt.201 zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)202 203 zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &204 & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )205 206 zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )207 208 ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &209 & + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &210 & +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)211 END DO212 END DO213 212 ENDIF 214 213 … … 216 215 DO jj = 1, jpjm1 217 216 DO ji = 1, fs_jpim1 ! vector opt. 218 zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) /e2f(ji,jj)219 220 zmskf = 1. /MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1)&221 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1.)222 223 zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )217 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * fse3f(ji,jj,jk) * r1_e2f(ji,jj) 218 219 zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & 220 & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ) , 1._wp ) 221 222 zcof2 = - rn_aht_0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) 224 223 225 224 zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) & … … 237 236 DO jj = 2, jpjm1 238 237 DO ji = 1, fs_jpim1 ! vector opt. 239 zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) /e1f(ji,jj)240 241 zmskf = 1. /MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1)&242 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1.)243 244 zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )245 246 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &247 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)&248 & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)238 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * fse3f(ji,jj,jk) * r1_e1f(ji,jj) 239 240 zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & 241 & + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 242 243 zcof1 = - rn_aht_0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) ) 244 245 zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) & 246 & + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) & 247 & + zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk) 249 248 END DO 250 249 END DO … … 254 253 DO jj = 2, jpj 255 254 DO ji = 1, fs_jpim1 ! vector opt. 256 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj) 255 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) * r1_e2t(ji,jj) 256 257 zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 258 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ) , 1._wp ) 259 260 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 261 262 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & 263 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) & 264 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk) 265 END DO 266 END DO 267 ELSE ! other coordinate system (zco or sco) : e3t 268 DO jj = 2, jpj 269 DO ji = 1, fs_jpim1 ! vector opt. 270 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * fse3t(ji,jj,jk) * r1_e2t(ji,jj) 257 271 258 272 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & 259 273 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) 260 274 261 zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )275 zcof2 = - rn_aht_0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) ) 262 276 263 277 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) & … … 266 280 END DO 267 281 END DO 268 ELSE ! other coordinate system (zco or sco) : e3t269 DO jj = 2, jpj270 DO ji = 1, fs_jpim1 ! vector opt.271 zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)272 273 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &274 & + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )275 276 zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )277 278 zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &279 & + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &280 & +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)281 END DO282 END DO283 282 ENDIF 284 283 … … 286 285 ! Second derivative (divergence) and add to the general trend 287 286 ! ----------------------------------------------------------- 288 289 287 DO jj = 2, jpjm1 290 DO ji = 2, jpim1 !! Question vectop possible??? !!bug 291 ! volume elements 292 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 293 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 294 ! horizontal component of isopycnal momentum diffusive trends 295 zuah =( ziut (ji+1,jj) - ziut (ji,jj ) + & 296 & zjuf (ji ,jj) - zjuf (ji,jj-1) ) / zbu 297 zvah =( zivf (ji,jj ) - zivf (ji-1,jj) + & 298 & zjvt (ji,jj+1) - zjvt (ji,jj ) ) / zbv 299 ! add the trends to the general trends 300 ua (ji,jj,jk) = ua (ji,jj,jk) + zuah 301 va (ji,jj,jk) = va (ji,jj,jk) + zvah 288 DO ji = 2, jpim1 !!gm Question vectop possible??? !!bug 289 ua(ji,jj,jk) = ua(ji,jj,jk) + ( ziut(ji+1,jj) - ziut(ji,jj ) & 290 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 291 va(ji,jj,jk) = va(ji,jj,jk) + ( zivf(ji,jj ) - zivf(ji-1,jj) & 292 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 302 293 END DO 303 294 END DO … … 358 349 DO jk = 2, jpkm1 359 350 DO ji = 2, jpim1 360 zcoef0= 0.5 * aht0 * umask(ji,jj,jk)361 351 zcoef0= 0.5 * rn_aht_0 * umask(ji,jj,jk) 352 ! 362 353 zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) ) 363 354 zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) ) 364 355 ! 365 356 zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) & 366 357 + tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. ) 367 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) +fmask(ji,jj,jk-1) &368 + fmask(ji,jj-1,jk ) +fmask(ji,jj,jk ), 1. )358 zmkf = 1./MAX( fmask(ji,jj-1,jk-1) + fmask(ji,jj,jk-1) & 359 + fmask(ji,jj-1,jk ) + fmask(ji,jj,jk ), 1. ) 369 360 370 361 zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi … … 376 367 +zdj1u(ji,jk ) + zdju (ji ,jk ) ) 377 368 ! update avmu (add isopycnal vertical coefficient to avmu) 378 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0379 avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0369 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 370 avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / rn_aht_0 380 371 END DO 381 372 END DO … … 384 375 DO jk = 2, jpkm1 385 376 DO ji = 2, jpim1 386 zcoef0 = 0.5 * aht0 * vmask(ji,jj,jk)377 zcoef0 = 0.5 * rn_aht_0 * vmask(ji,jj,jk) 387 378 388 379 zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) ) … … 398 389 ! vertical flux on v field 399 390 zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) & 400 401 402 391 & +zdiv (ji,jk ) + zdiv (ji-1,jk ) ) & 392 & + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) & 393 & +zdjv (ji,jk ) + zdj1v(ji ,jk ) ) 403 394 ! update avmv (add isopycnal vertical coefficient to avmv) 404 ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0405 avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0395 ! Caution: zcoef0 include rn_aht_0, so divided by rn_aht_0 to obtain slp^2 * rn_aht_0 396 avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / rn_aht_0 406 397 END DO 407 398 END DO … … 412 403 DO jk = 1, jpkm1 413 404 DO ji = 2, jpim1 414 ! volume elements 415 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 416 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 417 ! part of the k-component of isopycnal momentum diffusive trends 418 zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu 419 zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv 420 ! add the trends to the general trends 421 ua(ji,jj,jk) = ua(ji,jj,jk) + zuav 422 va(ji,jj,jk) = va(ji,jj,jk) + zvav 405 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 406 va(ji,jj,jk) = va(ji,jj,jk) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 423 407 END DO 424 408 END DO … … 432 416 END SUBROUTINE dyn_ldf_iso 433 417 434 # else435 !!----------------------------------------------------------------------436 !! Dummy module NO rotation of mixing tensor437 !!----------------------------------------------------------------------438 CONTAINS439 SUBROUTINE dyn_ldf_iso( kt ) ! Empty routine440 INTEGER, INTENT(in) :: kt441 WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt442 END SUBROUTINE dyn_ldf_iso443 #endif444 445 418 !!====================================================================== 446 419 END MODULE dynldf_iso -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5120 r5836 77 77 !! Note that as all external forcing a time averaging over a two rdt 78 78 !! period is used to prevent the divergence of odd and even time step. 79 !!80 !! N.B. : When key_esopa is used all the scheme are tested, regardless81 !! of the physical meaning of the results.82 79 !!---------------------------------------------------------------------- 83 80 ! … … 121 118 DO ji = fs_2, fs_jpim1 ! vector opt. 122 119 spgu(ji,jj) = spgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) & 123 & + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) /e1u(ji,jj)120 & + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 124 121 spgv(ji,jj) = spgv(ji,jj) + zg_2 * ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) & 125 & + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)122 & + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 126 123 END DO 127 124 END DO … … 135 132 DO jj = 2, jpjm1 ! add tide potential forcing 136 133 DO ji = fs_2, fs_jpim1 ! vector opt. 137 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) /e1u(ji,jj)138 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) /e2v(ji,jj)134 spgu(ji,jj) = spgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 135 spgv(ji,jj) = spgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 139 136 END DO 140 137 END DO … … 149 146 DO jj = 2, jpjm1 150 147 DO ji = fs_2, fs_jpim1 ! vector opt. 151 spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) /e1u(ji,jj)152 spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) /e2v(ji,jj)148 spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 149 spgv(ji,jj) = spgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 153 150 END DO 154 151 END DO … … 176 173 CASE ( 2 ) ; CALL dyn_spg_flt( kt, kindic ) ! filtered 177 174 ! 178 CASE ( -1 ) ! esopa: test all possibility with control print179 CALL dyn_spg_exp( kt )180 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg0 - Ua: ', mask1=umask, &181 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )182 CALL dyn_spg_ts ( kt )183 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg1 - Ua: ', mask1=umask, &184 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )185 CALL dyn_spg_flt( kt, kindic )186 CALL prt_ctl( tab3d_1=ua, clinfo1=' spg2 - Ua: ', mask1=umask, &187 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )188 175 END SELECT 189 176 ! … … 248 235 IF(lk_dynspg_flt) ioptio = ioptio + 1 249 236 ! 250 IF( ( ioptio > 1 .AND. .NOT. lk_esopa ).OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) &237 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 251 238 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 252 239 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. ln_isfcav ) & 253 240 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 254 241 ! 255 IF( lk_esopa ) nspg = -1256 242 IF( lk_dynspg_exp) nspg = 0 257 243 IF( lk_dynspg_ts ) nspg = 1 258 244 IF( lk_dynspg_flt) nspg = 2 259 245 ! 260 IF( lk_esopa ) nspg = -1261 !262 246 IF(lwp) THEN 263 247 WRITE(numout,*) 264 IF( nspg == -1 ) WRITE(numout,*) ' ESOPA test All scheme used'265 248 IF( nspg == 0 ) WRITE(numout,*) ' explicit free surface' 266 249 IF( nspg == 1 ) WRITE(numout,*) ' free surface with time splitting scheme' … … 268 251 ENDIF 269 252 270 #if defined key_dynspg_flt || defined key_esopa253 #if defined key_dynspg_flt 271 254 CALL solver_init( nit000 ) ! Elliptic solver initialisation 272 255 #endif 273 274 ! ! Control of timestep choice275 IF( lk_dynspg_ts .OR. lk_dynspg_exp ) THEN276 IF( nn_cla == 1 ) CALL ctl_stop( 'Crossland advection not implemented for this free surface formulation' )277 ENDIF278 279 256 ! ! Control of hydrostatic pressure choice 280 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN 281 CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 282 ENDIF 257 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 283 258 ! 284 259 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_init') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r4990 r5836 7 7 !! 3.2 ! 2009-06 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 8 8 !!---------------------------------------------------------------------- 9 #if defined key_dynspg_exp || defined key_esopa9 #if defined key_dynspg_exp 10 10 !!---------------------------------------------------------------------- 11 11 !! 'key_dynspg_exp' explicit free surface -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r4990 r5836 14 14 !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 15 15 !! 3.7 ! 2014-04 (F. Roquet, G. Madec) add some trends diag 16 !!---------------------------------------------------------------------- 17 #if defined key_dynspg_flt || defined key_esopa 16 !! - ! 2014-12 (G. Madec) remove cross-land advection (cla) 17 !!---------------------------------------------------------------------- 18 #if defined key_dynspg_flt 18 19 !!---------------------------------------------------------------------- 19 20 !! 'key_dynspg_flt' filtered free surface … … 36 37 USE bdydyn ! ocean open boundary condition on dynamics 37 38 USE bdyvol ! ocean open boundary condition (bdy_vol routine) 38 USE cla ! cross land advection39 39 USE trd_oce ! trends: ocean variables 40 40 USE trddyn ! trend manager: dynamics … … 206 206 CALL Agrif_dyn( kt ) ! Update velocities on each coarse/fine interfaces 207 207 #endif 208 IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 ) CALL cla_dynspg( kt ) ! Cross Land Advection (update (ua,va))209 208 210 209 ! compute the next vertically averaged velocity (effect of the additional force not included) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_oce.F90
r4486 r5836 17 17 18 18 ! !!! Surface pressure gradient logicals 19 #if defined key_dynspg_exp || defined key_esopa19 #if defined key_dynspg_exp 20 20 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_exp = .TRUE. !: Explicit free surface flag 21 21 #else 22 22 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_exp = .FALSE. !: Explicit free surface flag 23 23 #endif 24 #if defined key_dynspg_ts || defined key_esopa24 #if defined key_dynspg_ts 25 25 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_ts = .TRUE. !: Free surface with time splitting flag 26 26 #else 27 27 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_ts = .FALSE. !: Free surface with time splitting flag 28 28 #endif 29 #if defined key_dynspg_flt || defined key_esopa29 #if defined key_dynspg_flt 30 30 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_flt = .TRUE. !: Filtered free surface cst volume flag 31 31 #else 32 32 LOGICAL, PUBLIC, PARAMETER :: lk_dynspg_flt = .FALSE. !: Filtered free surface cst volume flag 33 33 #endif 34 35 34 ! !!! Time splitting scheme (key_dynspg_ts) 36 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshn_e, ssha_e ! sea surface heigth (now, after) -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5656 r5836 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 13 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts || defined key_esopa14 #if defined key_dynspg_ts 15 15 !!---------------------------------------------------------------------- 16 16 !! 'key_dynspg_ts' split explicit free surface … … 98 98 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 99 99 100 IF( ln_dynvor_een .or. ln_dynvor_een_old) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , &101 &ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) )100 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 101 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 102 102 103 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) … … 107 107 ! 108 108 END FUNCTION dyn_spg_ts_alloc 109 109 110 110 111 SUBROUTINE dyn_spg_ts( kt ) … … 219 220 ! 220 221 IF ( kt == nit000 .OR. lk_vvl ) THEN 221 IF ( ln_dynvor_een_old ) THEN 222 DO jj = 1, jpjm1 223 DO ji = 1, jpim1 224 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 225 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 226 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 227 END DO 228 END DO 222 IF ( ln_dynvor_een ) THEN !== EEN scheme ==! 223 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 224 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 225 DO jj = 1, jpjm1 226 DO ji = 1, jpim1 227 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 228 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 229 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 230 END DO 231 END DO 232 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 233 DO jj = 1, jpjm1 234 DO ji = 1, jpim1 235 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 236 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 237 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 238 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 239 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 240 END DO 241 END DO 242 END SELECT 229 243 CALL lbc_lnk( zwz, 'F', 1._wp ) 230 zwz(:,:) = ff(:,:) * zwz(:,:) 231 244 ! 232 245 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 233 246 DO jj = 2, jpj 234 DO ji = fs_2, jpi ! vector opt.247 DO ji = 2, jpi 235 248 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 236 249 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 239 252 END DO 240 253 END DO 241 ELSE IF ( ln_dynvor_een ) THEN 242 DO jj = 1, jpjm1 243 DO ji = 1, jpim1 244 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 245 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 246 & / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 247 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 248 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 249 END DO 250 END DO 251 CALL lbc_lnk( zwz, 'F', 1._wp ) 252 zwz(:,:) = ff(:,:) * zwz(:,:) 253 254 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 255 DO jj = 2, jpj 256 DO ji = fs_2, jpi ! vector opt. 257 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 258 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 259 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 260 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 261 END DO 262 END DO 263 ELSE 254 ! 255 ELSE !== all other schemes (ENE, ENS, MIX) 264 256 zwz(:,:) = 0._wp 265 zhf(:,:) = 0. 257 zhf(:,:) = 0._wp 266 258 IF ( .not. ln_sco ) THEN 259 260 !!gm agree the JC comment : this should be done in a much clear way 261 267 262 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 268 263 ! Set it to zero for the time being … … 276 271 277 272 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) *(1._wp- umask(:,jj,1) * umask(:,jj+1,1))273 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 279 274 END DO 280 275 … … 297 292 ! If forward start at previous time step, and centered integration, 298 293 ! then update averaging weights: 299 IF ( (.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN294 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 300 295 ll_fw_start=.FALSE. 301 296 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) … … 338 333 DO jj = 2, jpjm1 339 334 DO ji = fs_2, fs_jpim1 ! vector opt. 340 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)341 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)342 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) /e2v(ji,jj)343 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) /e2v(ji,jj)335 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 336 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 337 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 338 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 344 339 ! energy conserving formulation for planetary vorticity term 345 340 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) … … 352 347 DO ji = fs_2, fs_jpim1 ! vector opt. 353 348 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 354 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)349 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 355 350 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 356 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)351 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 357 352 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 358 353 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 360 355 END DO 361 356 ! 362 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN ! enstrophy and energy conserving scheme357 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 363 358 DO jj = 2, jpjm1 364 359 DO ji = fs_2, fs_jpim1 ! vector opt. 365 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &366 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &367 & + ftse(ji,jj ) * zwy(ji ,jj-1) &368 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )369 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &370 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &371 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &372 & + ftne(ji,jj ) * zwx(ji ,jj ) )360 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 361 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 362 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 363 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 364 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 365 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 366 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 367 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 373 368 END DO 374 369 END DO … … 381 376 DO jj = 2, jpjm1 382 377 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) /e1u(ji,jj)384 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) /e2v(ji,jj)378 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 379 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 385 380 END DO 386 381 END DO … … 431 426 DO jj = 2, jpjm1 432 427 DO ji = fs_2, fs_jpim1 ! vector opt. 433 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)434 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)428 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 429 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 435 430 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 436 431 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 441 436 DO ji = fs_2, fs_jpim1 ! vector opt. 442 437 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 443 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)438 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 444 439 zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 445 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)440 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 446 441 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 447 442 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 549 544 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 550 545 DO ji = 2, fs_jpim1 ! Vector opt. 551 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1 2u(ji,jj) &552 & * ( e1 2t(ji ,jj) * zsshp2_e(ji ,jj) &553 & + e1 2t(ji+1,jj) * zsshp2_e(ji+1,jj) )554 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1 2v(ji,jj) &555 & * ( e1 2t(ji,jj ) * zsshp2_e(ji,jj ) &556 & + e1 2t(ji,jj+1) * zsshp2_e(ji,jj+1) )546 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 547 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 548 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 549 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 550 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 551 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 557 552 END DO 558 553 END DO … … 602 597 ! Sum over sub-time-steps to compute advective velocities 603 598 za2 = wgtbtp2(jn) 604 zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u(:,:)605 zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v(:,:)599 zu_sum(:,:) = zu_sum(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 600 zv_sum(:,:) = zv_sum(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 606 601 ! 607 602 ! Set next sea level: … … 609 604 DO ji = fs_2, fs_jpim1 ! vector opt. 610 605 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 611 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1 2t(ji,jj)606 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 612 607 END DO 613 608 END DO … … 627 622 DO jj = 2, jpjm1 628 623 DO ji = 2, jpim1 ! NO Vector Opt. 629 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1 2u(ji,jj) &630 & * ( e1 2t(ji ,jj ) * ssha_e(ji ,jj ) &631 & + e1 2t(ji+1,jj ) * ssha_e(ji+1,jj ) )632 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1 2v(ji,jj) &633 & * ( e1 2t(ji ,jj ) * ssha_e(ji ,jj ) &634 & + e1 2t(ji ,jj+1) * ssha_e(ji ,jj+1) )624 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 625 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 626 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 627 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 628 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 629 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 635 630 END DO 636 631 END DO … … 666 661 DO jj = 2, jpjm1 667 662 DO ji = 2, jpim1 668 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1 2u(ji ,jj) &669 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj) &670 & + e1 2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )671 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1 2v(ji ,jj ) &672 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj ) &673 & + e1 2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )663 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1e2u(ji ,jj) & 664 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 665 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 666 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1e2v(ji ,jj ) & 667 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 668 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 674 669 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 675 670 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 … … 688 683 DO jj = 2, jpjm1 689 684 DO ji = fs_2, fs_jpim1 ! vector opt. 690 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)691 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)692 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) /e2v(ji,jj)693 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)685 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 686 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 687 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 688 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 694 689 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 695 690 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) … … 701 696 DO ji = fs_2, fs_jpim1 ! vector opt. 702 697 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 703 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)698 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 704 699 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 705 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)700 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 706 701 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 707 702 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 709 704 END DO 710 705 ! 711 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN !== energy and enstrophy conserving scheme ==!706 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 712 707 DO jj = 2, jpjm1 713 708 DO ji = fs_2, fs_jpim1 ! vector opt. 714 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &715 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &716 & + ftse(ji,jj ) * zwy(ji ,jj-1) &717 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )718 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &719 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &720 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &721 & + ftne(ji,jj ) * zwx(ji ,jj ) )709 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 710 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 711 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 712 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 713 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 714 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 715 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 716 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 722 717 END DO 723 718 END DO … … 729 724 DO jj = 2, jpjm1 730 725 DO ji = fs_2, fs_jpim1 ! vector opt. 731 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) /e1u(ji,jj)732 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) /e2v(ji,jj)726 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 727 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 733 728 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 734 729 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg … … 745 740 DO ji = fs_2, fs_jpim1 ! vector opt. 746 741 ! Add surface pressure gradient 747 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) /e1u(ji,jj)748 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) /e2v(ji,jj)742 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 743 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 749 744 zwx(ji,jj) = zu_spg 750 745 zwy(ji,jj) = zv_spg … … 850 845 DO jj = 1, jpjm1 851 846 DO ji = 1, jpim1 ! NO Vector Opt. 852 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1 2u(ji,jj) &853 & * ( e1 2t(ji ,jj) * ssha(ji ,jj) &854 & + e1 2t(ji+1,jj) * ssha(ji+1,jj) )855 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1 2v(ji,jj) &856 & * ( e1 2t(ji,jj ) * ssha(ji,jj ) &857 & + e1 2t(ji,jj+1) * ssha(ji,jj+1) )847 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 848 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 849 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 850 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 851 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 852 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 858 853 END DO 859 854 END DO … … 1093 1088 DO jj = 1, jpj 1094 1089 DO ji =1, jpi 1095 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))1096 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))1097 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) )1090 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1091 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1092 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1098 1093 END DO 1099 1094 END DO … … 1101 1096 DO jj = 1, jpj 1102 1097 DO ji =1, jpi 1103 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj))1104 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj))1105 zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) )1106 END DO 1107 END DO 1108 ENDIF 1109 1110 zcmax = MAXVAL( zcu(:,:))1098 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1099 zyr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1100 zcu(ji,jj) = SQRT( grav * ht(ji,jj) * (zxr2 + zyr2) ) 1101 END DO 1102 END DO 1103 ENDIF 1104 1105 zcmax = MAXVAL( zcu(:,:) ) 1111 1106 IF( lk_mpp ) CALL mpp_max( zcmax ) 1112 1107 … … 1114 1109 IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1115 1110 1116 rdtbt = rdt / FLOAT(nn_baro)1111 rdtbt = rdt / REAL( nn_baro , wp ) 1117 1112 zcmax = zcmax * rdtbt 1118 1113 ! Print results … … 1194 1189 !!====================================================================== 1195 1190 END MODULE dynspg_ts 1196 1197 1198 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r5029 r5836 16 16 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 17 17 !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 18 !! - ! 2014-06 (G. Madec) suppression of velocity curl from in-core memory 18 19 !!---------------------------------------------------------------------- 19 20 … … 22 23 !! vor_ens : enstrophy conserving scheme (ln_dynvor_ens=T) 23 24 !! vor_ene : energy conserving scheme (ln_dynvor_ene=T) 24 !! vor_mix : mixed enstrophy/energy conserving (ln_dynvor_mix=T)25 25 !! vor_een : energy and enstrophy conserving (ln_dynvor_een=T) 26 26 !! dyn_vor_init : set and control of the different vorticity option … … 32 32 USE trd_oce ! trends: ocean variables 33 33 USE trddyn ! trend manager: dynamics 34 ! 34 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 36 USE prtctl ! Print control … … 44 45 45 46 PUBLIC dyn_vor ! routine called by step.F90 46 PUBLIC dyn_vor_init ! routine called by opa.F9047 PUBLIC dyn_vor_init ! routine called by nemogcm.F90 47 48 48 49 ! !!* Namelist namdyn_vor: vorticity term 49 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme 50 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme 51 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme 52 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme 53 LOGICAL, PUBLIC :: ln_dynvor_een_old !: energy and enstrophy conserving scheme (original formulation) 54 55 INTEGER :: nvor = 0 ! type of vorticity trend used 56 INTEGER :: ncor = 1 ! coriolis 57 INTEGER :: nrvm = 2 ! =2 relative vorticity ; =3 metric term 58 INTEGER :: ntot = 4 ! =4 total vorticity (relative + planetary) ; =5 coriolis + metric term 59 50 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme (ENE) 51 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 52 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 53 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme (EEN) 54 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 55 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 56 57 INTEGER :: nvor_scheme ! choice of the type of advection scheme 58 ! ! associated indices: 59 INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme 60 INTEGER, PUBLIC, PARAMETER :: np_ENS = 2 ! ENS scheme 61 INTEGER, PUBLIC, PARAMETER :: np_MIX = 3 ! MIX scheme 62 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 63 64 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 65 ! ! associated indices: 66 INTEGER, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 67 INTEGER, PARAMETER :: np_RVO = 2 ! relative vorticity 68 INTEGER, PARAMETER :: np_MET = 3 ! metric term 69 INTEGER, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 70 INTEGER, PARAMETER :: np_CME = 5 ! Coriolis + metric term 71 72 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 73 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 74 REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12 75 60 76 !! * Substitutions 61 77 # include "domzgr_substitute.h90" 62 78 # include "vectopt_loop_substitute.h90" 63 79 !!---------------------------------------------------------------------- 64 !! NEMO/OPA 3. 3 , NEMO Consortium (2010)80 !! NEMO/OPA 3.7 , NEMO Consortium (2014) 65 81 !! $Id$ 66 82 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 87 103 IF( l_trddyn ) CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 88 104 ! 89 ! ! vorticity term 90 SELECT CASE ( nvor ) ! compute the vorticity trend and add it to the general trend 91 ! 92 CASE ( -1 ) ! esopa: test all possibility with control print 93 CALL vor_ene( kt, ntot, ua, va ) 94 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor0 - Ua: ', mask1=umask, & 95 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 96 CALL vor_ens( kt, ntot, ua, va ) 97 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor1 - Ua: ', mask1=umask, & 98 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 99 CALL vor_mix( kt ) 100 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor2 - Ua: ', mask1=umask, & 101 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 102 CALL vor_een( kt, ntot, ua, va ) 103 CALL prt_ctl( tab3d_1=ua, clinfo1=' vor3 - Ua: ', mask1=umask, & 104 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 105 ! 106 CASE ( 0 ) ! energy conserving scheme 107 IF( l_trddyn ) THEN 108 ztrdu(:,:,:) = ua(:,:,:) 109 ztrdv(:,:,:) = va(:,:,:) 110 CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend 105 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 106 ! 107 CASE ( np_ENE ) !* energy conserving scheme 108 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 109 ztrdu(:,:,:) = ua(:,:,:) 110 ztrdv(:,:,:) = va(:,:,:) 111 CALL vor_ene( kt, nrvm, ua, va ) ! relative vorticity or metric trend 111 112 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 112 113 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 114 115 ztrdu(:,:,:) = ua(:,:,:) 115 116 ztrdv(:,:,:) = va(:,:,:) 116 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend117 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend 117 118 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 118 119 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 119 120 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 120 121 ELSE 121 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity 122 ENDIF 123 ! 124 CASE ( 1 ) !enstrophy conserving scheme125 IF( l_trddyn ) THEN126 ztrdu(:,:,:) = ua(:,:,:) 127 ztrdv(:,:,:) = va(:,:,:) 128 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend122 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity trend 123 ENDIF 124 ! 125 CASE ( np_ENS ) !* enstrophy conserving scheme 126 IF( l_trddyn ) THEN ! trend diagnostics: splitthe trend in two 127 ztrdu(:,:,:) = ua(:,:,:) 128 ztrdv(:,:,:) = va(:,:,:) 129 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend 129 130 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 130 131 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 132 133 ztrdu(:,:,:) = ua(:,:,:) 133 134 ztrdv(:,:,:) = va(:,:,:) 134 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend135 CALL vor_ens( kt, ncor, ua, va ) ! planetary vorticity trend 135 136 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 136 137 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 137 138 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 138 139 ELSE 139 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity 140 ENDIF 141 ! 142 CASE ( 2 ) !mixed ene-ens scheme143 IF( l_trddyn ) THEN144 ztrdu(:,:,:) = ua(:,:,:) 145 ztrdv(:,:,:) = va(:,:,:) 146 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens)140 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity trend 141 ENDIF 142 ! 143 CASE ( np_MIX ) !* mixed ene-ens scheme 144 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 145 ztrdu(:,:,:) = ua(:,:,:) 146 ztrdv(:,:,:) = va(:,:,:) 147 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 147 148 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 148 149 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 150 151 ztrdu(:,:,:) = ua(:,:,:) 151 152 ztrdv(:,:,:) = va(:,:,:) 152 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene)153 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 153 154 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 154 155 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 155 156 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 157 ELSE 157 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) 158 ENDIF 159 ! 160 CASE ( 3 ) ! energy and enstrophy conserving scheme 161 IF( l_trddyn ) THEN 162 ztrdu(:,:,:) = ua(:,:,:) 163 ztrdv(:,:,:) = va(:,:,:) 164 CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend 158 CALL vor_ens( kt, nrvm, ua, va ) ! relative vorticity or metric trend (ens) 159 CALL vor_ene( kt, ncor, ua, va ) ! planetary vorticity trend (ene) 160 ENDIF 161 ! 162 CASE ( np_EEN ) !* energy and enstrophy conserving scheme 163 IF( l_trddyn ) THEN ! trend diagnostics: split the trend in two 164 ztrdu(:,:,:) = ua(:,:,:) 165 ztrdv(:,:,:) = va(:,:,:) 166 CALL vor_een( kt, nrvm, ua, va ) ! relative vorticity or metric trend 165 167 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 166 168 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) … … 168 170 ztrdu(:,:,:) = ua(:,:,:) 169 171 ztrdv(:,:,:) = va(:,:,:) 170 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend172 CALL vor_een( kt, ncor, ua, va ) ! planetary vorticity trend 171 173 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 174 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 175 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 174 176 ELSE 175 CALL vor_een( kt, ntot, ua, va ) ! total vorticity 177 CALL vor_een( kt, ntot, ua, va ) ! total vorticity trend 176 178 ENDIF 177 179 ! … … 197 199 !! 198 200 !! ** Method : Trend evaluated using now fields (centered in time) 199 !! and the Sadourny (1975) flux form formulation : conserves the 200 !! horizontal kinetic energy. 201 !! The trend of the vorticity term is given by: 202 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives: 203 !! voru = 1/e1u mj-1[ (rotn+f)/e3f mi(e1v*e3v vn) ] 204 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f mj(e2u*e3u un) ] 205 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 206 !! voru = 1/e1u mj-1[ (rotn+f) mi(e1v vn) ] 207 !! vorv = 1/e2v mi-1[ (rotn+f) mj(e2u un) ] 208 !! Add this trend to the general momentum trend (ua,va): 209 !! (ua,va) = (ua,va) + ( voru , vorv ) 201 !! and the Sadourny (1975) flux form formulation : conserves the 202 !! horizontal kinetic energy. 203 !! The general trend of momentum is increased due to the vorticity 204 !! term which is given by: 205 !! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v vn) ] 206 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f mj(e2u*e3u un) ] 207 !! where rvor is the relative vorticity 210 208 !! 211 209 !! ** Action : - Update (ua,va) with the now vorticity term trend … … 213 211 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 214 212 !!---------------------------------------------------------------------- 215 !216 213 INTEGER , INTENT(in ) :: kt ! ocean time-step index 217 214 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 220 217 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 221 218 ! 222 INTEGER :: ji, jj, jk ! dummy loop indices223 REAL(wp) :: zx1, zy1, z fact2, zx2, zy2 ! local scalars224 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz219 INTEGER :: ji, jj, jk ! dummy loop indices 220 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 221 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz ! 2D workspace 225 222 !!---------------------------------------------------------------------- 226 223 ! … … 234 231 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 235 232 ENDIF 236 237 zfact2 = 0.5 * 0.5 ! Local constant initialization 238 239 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz ) 233 ! 240 234 ! ! =============== 241 235 DO jk = 1, jpkm1 ! Horizontal slab 242 236 ! ! =============== 243 237 ! 244 ! Potential vorticity and horizontal fluxes 245 ! ----------------------------------------- 246 SELECT CASE( kvor ) ! vorticity considered 247 CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) 248 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity 249 CASE ( 3 ) ! metric term 238 SELECT CASE( kvor ) !== vorticity considered ==! 239 CASE ( np_COR ) !* Coriolis (planetary vorticity) 240 zwz(:,:) = ff(:,:) 241 CASE ( np_RVO ) !* relative vorticity 242 DO jj = 1, jpjm1 243 DO ji = 1, fs_jpim1 ! vector opt. 244 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 245 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 246 END DO 247 END DO 248 CASE ( np_MET ) !* metric term 250 249 DO jj = 1, jpjm1 251 250 DO ji = 1, fs_jpim1 ! vector opt. 252 251 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 253 252 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 254 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 255 END DO 256 END DO 257 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) 258 CASE ( 5 ) ! total (coriolis + metric) 259 DO jj = 1, jpjm1 260 DO ji = 1, fs_jpim1 ! vector opt. 261 zwz(ji,jj) = ( ff (ji,jj) & 262 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 263 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 264 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 265 & ) 266 END DO 267 END DO 253 & * 0.5 * r1_e1e2f(ji,jj) 254 END DO 255 END DO 256 CASE ( np_CRV ) !* Coriolis + relative vorticity 257 DO jj = 1, jpjm1 258 DO ji = 1, fs_jpim1 ! vector opt. 259 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 260 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 261 & * r1_e1e2f(ji,jj) 262 END DO 263 END DO 264 CASE ( np_CME ) !* Coriolis + metric 265 DO jj = 1, jpjm1 266 DO ji = 1, fs_jpim1 ! vector opt. 267 zwz(ji,jj) = ff(ji,jj) & 268 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 269 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 270 & * 0.5 * r1_e1e2f(ji,jj) 271 END DO 272 END DO 273 CASE DEFAULT ! error 274 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 268 275 END SELECT 276 ! 277 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 278 DO jj = 1, jpjm1 279 DO ji = 1, fs_jpim1 ! vector opt. 280 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 281 END DO 282 END DO 283 ENDIF 269 284 270 285 IF( ln_sco ) THEN … … 276 291 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 277 292 ENDIF 278 279 ! Compute and add the vorticity term trend 280 ! ---------------------------------------- 293 ! !== compute and add the vorticity term trend =! 281 294 DO jj = 2, jpjm1 282 295 DO ji = fs_2, fs_jpim1 ! vector opt. … … 285 298 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 286 299 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 287 pua(ji,jj,jk) = pua(ji,jj,jk) + zfact2 /e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )288 pva(ji,jj,jk) = pva(ji,jj,jk) - zfact2 /e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )300 pua(ji,jj,jk) = pua(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 301 pva(ji,jj,jk) = pva(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 289 302 END DO 290 303 END DO … … 299 312 300 313 301 SUBROUTINE vor_mix( kt )302 !!----------------------------------------------------------------------303 !! *** ROUTINE vor_mix ***304 !!305 !! ** Purpose : Compute the now total vorticity trend and add it to306 !! the general trend of the momentum equation.307 !!308 !! ** Method : Trend evaluated using now fields (centered in time)309 !! Mixte formulation : conserves the potential enstrophy of a hori-310 !! zontally non-divergent flow for (rotzu x uh), the relative vor-311 !! ticity term and the horizontal kinetic energy for (f x uh), the312 !! coriolis term. the now trend of the vorticity term is given by:313 !! * s-coordinate (ln_sco=T), the e3. are inside the derivatives:314 !! voru = 1/e1u mj-1(rotn/e3f) mj-1[ mi(e1v*e3v vn) ]315 !! +1/e1u mj-1[ f/e3f mi(e1v*e3v vn) ]316 !! vorv = 1/e2v mi-1(rotn/e3f) mi-1[ mj(e2u*e3u un) ]317 !! +1/e2v mi-1[ f/e3f mj(e2u*e3u un) ]318 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes:319 !! voru = 1/e1u mj-1(rotn) mj-1[ mi(e1v vn) ]320 !! +1/e1u mj-1[ f mi(e1v vn) ]321 !! vorv = 1/e2v mi-1(rotn) mi-1[ mj(e2u un) ]322 !! +1/e2v mi-1[ f mj(e2u un) ]323 !! Add this now trend to the general momentum trend (ua,va):324 !! (ua,va) = (ua,va) + ( voru , vorv )325 !!326 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend327 !!328 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689.329 !!----------------------------------------------------------------------330 !331 INTEGER, INTENT(in) :: kt ! ocean timestep index332 !333 INTEGER :: ji, jj, jk ! dummy loop indices334 REAL(wp) :: zfact1, zua, zcua, zx1, zy1 ! local scalars335 REAL(wp) :: zfact2, zva, zcva, zx2, zy2 ! - -336 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww337 !!----------------------------------------------------------------------338 !339 IF( nn_timing == 1 ) CALL timing_start('vor_mix')340 !341 CALL wrk_alloc( jpi, jpj, zwx, zwy, zwz, zww )342 !343 IF( kt == nit000 ) THEN344 IF(lwp) WRITE(numout,*)345 IF(lwp) WRITE(numout,*) 'dyn:vor_mix : vorticity term: mixed energy/enstrophy conserving scheme'346 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'347 ENDIF348 349 zfact1 = 0.5 * 0.25 ! Local constant initialization350 zfact2 = 0.5 * 0.5351 352 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, zww )353 ! ! ===============354 DO jk = 1, jpkm1 ! Horizontal slab355 ! ! ===============356 !357 ! Relative and planetary potential vorticity and horizontal fluxes358 ! ----------------------------------------------------------------359 IF( ln_sco ) THEN360 IF( ln_dynadv_vec ) THEN361 zww(:,:) = rotn(:,:,jk) / fse3f(:,:,jk)362 ELSE363 DO jj = 1, jpjm1364 DO ji = 1, fs_jpim1 ! vector opt.365 zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &366 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &367 & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) * fse3f(ji,jj,jk) )368 END DO369 END DO370 ENDIF371 zwz(:,:) = ff (:,:) / fse3f(:,:,jk)372 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)373 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)374 ELSE375 IF( ln_dynadv_vec ) THEN376 zww(:,:) = rotn(:,:,jk)377 ELSE378 DO jj = 1, jpjm1379 DO ji = 1, fs_jpim1 ! vector opt.380 zww(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) &381 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) &382 & * 0.5 / ( e1f(ji,jj) * e2f (ji,jj) )383 END DO384 END DO385 ENDIF386 zwz(:,:) = ff (:,:)387 zwx(:,:) = e2u(:,:) * un(:,:,jk)388 zwy(:,:) = e1v(:,:) * vn(:,:,jk)389 ENDIF390 391 ! Compute and add the vorticity term trend392 ! ----------------------------------------393 DO jj = 2, jpjm1394 DO ji = fs_2, fs_jpim1 ! vector opt.395 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)396 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)397 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)398 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)399 ! enstrophy conserving formulation for relative vorticity term400 zua = zfact1 * ( zww(ji ,jj-1) + zww(ji,jj) ) * ( zy1 + zy2 )401 zva =-zfact1 * ( zww(ji-1,jj ) + zww(ji,jj) ) * ( zx1 + zx2 )402 ! energy conserving formulation for planetary vorticity term403 zcua = zfact2 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 )404 zcva =-zfact2 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 )405 ! mixed vorticity trend added to the momentum trends406 ua(ji,jj,jk) = ua(ji,jj,jk) + zcua + zua407 va(ji,jj,jk) = va(ji,jj,jk) + zcva + zva408 END DO409 END DO410 ! ! ===============411 END DO ! End of slab412 ! ! ===============413 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zwz, zww )414 !415 IF( nn_timing == 1 ) CALL timing_stop('vor_mix')416 !417 END SUBROUTINE vor_mix418 419 420 314 SUBROUTINE vor_ens( kt, kvor, pua, pva ) 421 315 !!---------------------------------------------------------------------- … … 429 323 !! potential enstrophy of a horizontally non-divergent flow. the 430 324 !! trend of the vorticity term is given by: 431 !! * s-coordinate (ln_sco=T), the e3. are inside the derivative: 432 !! voru = 1/e1u mj-1[ (rotn+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 433 !! vorv = 1/e2v mi-1[ (rotn+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 434 !! * z-coordinate (default key), e3t=e3u=e3v, the trend becomes: 435 !! voru = 1/e1u mj-1[ rotn+f ] mj-1[ mi(e1v vn) ] 436 !! vorv = 1/e2v mi-1[ rotn+f ] mi-1[ mj(e2u un) ] 325 !! voru = 1/e1u mj-1[ (rvor+f)/e3f ] mj-1[ mi(e1v*e3v vn) ] 326 !! vorv = 1/e2v mi-1[ (rvor+f)/e3f ] mi-1[ mj(e2u*e3u un) ] 437 327 !! Add this trend to the general momentum trend (ua,va): 438 328 !! (ua,va) = (ua,va) + ( voru , vorv ) … … 442 332 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 443 333 !!---------------------------------------------------------------------- 444 !445 334 INTEGER , INTENT(in ) :: kt ! ocean time-step index 446 335 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; … … 449 338 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 450 339 ! 451 INTEGER :: ji, jj, jk 452 REAL(wp) :: z fact1, zuav, zvau ! temporaryscalars453 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww340 INTEGER :: ji, jj, jk ! dummy loop indices 341 REAL(wp) :: zuav, zvau ! local scalars 342 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, zww ! 2D workspace 454 343 !!---------------------------------------------------------------------- 455 344 ! … … 463 352 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 464 353 ENDIF 465 466 zfact1 = 0.5 * 0.25 ! Local constant initialization467 468 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz )469 354 ! ! =============== 470 355 DO jk = 1, jpkm1 ! Horizontal slab 471 356 ! ! =============== 472 357 ! 473 ! Potential vorticity and horizontal fluxes 474 ! ----------------------------------------- 475 SELECT CASE( kvor ) ! vorticity considered 476 CASE ( 1 ) ; zwz(:,:) = ff(:,:) ! planetary vorticity (Coriolis) 477 CASE ( 2 ) ; zwz(:,:) = rotn(:,:,jk) ! relative vorticity 478 CASE ( 3 ) ! metric term 358 SELECT CASE( kvor ) !== vorticity considered ==! 359 CASE ( np_COR ) !* Coriolis (planetary vorticity) 360 zwz(:,:) = ff(:,:) 361 CASE ( np_RVO ) !* relative vorticity 362 DO jj = 1, jpjm1 363 DO ji = 1, fs_jpim1 ! vector opt. 364 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 365 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) * r1_e1e2f(ji,jj) 366 END DO 367 END DO 368 CASE ( np_MET ) !* metric term 479 369 DO jj = 1, jpjm1 480 370 DO ji = 1, fs_jpim1 ! vector opt. 481 371 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 482 372 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 483 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) 484 END DO 485 END DO 486 CASE ( 4 ) ; zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) ! total (relative + planetary vorticity) 487 CASE ( 5 ) ! total (coriolis + metric) 488 DO jj = 1, jpjm1 489 DO ji = 1, fs_jpim1 ! vector opt. 490 zwz(ji,jj) = ( ff (ji,jj) & 491 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 492 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 493 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 494 & ) 495 END DO 496 END DO 373 & * 0.5 * r1_e1e2f(ji,jj) 374 END DO 375 END DO 376 CASE ( np_CRV ) !* Coriolis + relative vorticity 377 DO jj = 1, jpjm1 378 DO ji = 1, fs_jpim1 ! vector opt. 379 zwz(ji,jj) = ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 380 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 381 & * r1_e1e2f(ji,jj) 382 END DO 383 END DO 384 CASE ( np_CME ) !* Coriolis + metric 385 DO jj = 1, jpjm1 386 DO ji = 1, fs_jpim1 ! vector opt. 387 zwz(ji,jj) = ff(ji,jj) & 388 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 389 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 390 & * 0.5 * r1_e1e2f(ji,jj) 391 END DO 392 END DO 393 CASE DEFAULT ! error 394 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 497 395 END SELECT 498 396 ! 499 IF( ln_sco ) THEN 500 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 501 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 502 zwz(ji,jj) = zwz(ji,jj) / fse3f(ji,jj,jk) 503 zwx(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk) * un(ji,jj,jk) 504 zwy(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk) * vn(ji,jj,jk) 505 END DO 506 END DO 397 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 398 DO jj = 1, jpjm1 399 DO ji = 1, fs_jpim1 ! vector opt. 400 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 401 END DO 402 END DO 403 ENDIF 404 ! 405 IF( ln_sco ) THEN !== horizontal fluxes ==! 406 zwz(:,:) = zwz(:,:) / fse3f(:,:,jk) 407 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 408 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 507 409 ELSE 508 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 509 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 510 zwx(ji,jj) = e2u(ji,jj) * un(ji,jj,jk) 511 zwy(ji,jj) = e1v(ji,jj) * vn(ji,jj,jk) 512 END DO 513 END DO 514 ENDIF 515 ! 516 ! Compute and add the vorticity term trend 517 ! ---------------------------------------- 410 zwx(:,:) = e2u(:,:) * un(:,:,jk) 411 zwy(:,:) = e1v(:,:) * vn(:,:,jk) 412 ENDIF 413 ! !== compute and add the vorticity term trend =! 518 414 DO jj = 2, jpjm1 519 415 DO ji = fs_2, fs_jpim1 ! vector opt. 520 zuav = zfact1 /e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) &521 & 522 zvau =- zfact1 /e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) &523 & 416 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 417 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 418 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 419 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 524 420 pua(ji,jj,jk) = pua(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 525 421 pva(ji,jj,jk) = pva(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 553 449 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 554 450 !!---------------------------------------------------------------------- 555 !556 451 INTEGER , INTENT(in ) :: kt ! ocean time-step index 557 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; 558 ! ! =nrvm (relative vorticity or metric) 452 INTEGER , INTENT(in ) :: kvor ! =ncor (planetary) ; =ntot (total) ; =nrvm (relative or metric) 559 453 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pua ! total u-trend 560 454 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pva ! total v-trend 561 !! 562 INTEGER :: ji, jj, jk ! dummy loop indices 563 INTEGER :: ierr ! local integer 564 REAL(wp) :: zfac12, zua, zva ! local scalars 565 REAL(wp) :: zmsk, ze3 ! local scalars 566 ! ! 3D workspace 567 REAL(wp), POINTER , DIMENSION(:,: ) :: zwx, zwy, zwz 568 REAL(wp), POINTER , DIMENSION(:,: ) :: ztnw, ztne, ztsw, ztse 569 #if defined key_vvl 570 REAL(wp), POINTER , DIMENSION(:,:,:) :: ze3f ! 3D workspace (lk_vvl=T) 571 #else 572 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:), SAVE :: ze3f ! lk_vvl=F, ze3f=1/e3f saved one for all 573 #endif 455 ! 456 INTEGER :: ji, jj, jk ! dummy loop indices 457 INTEGER :: ierr ! local integer 458 REAL(wp) :: zua, zva ! local scalars 459 REAL(wp) :: zmsk, ze3 ! local scalars 460 ! 461 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zwz, z1_e3f 462 REAL(wp), POINTER, DIMENSION(:,:) :: ztnw, ztne, ztsw, ztse 574 463 !!---------------------------------------------------------------------- 575 464 ! 576 465 IF( nn_timing == 1 ) CALL timing_start('vor_een') 577 466 ! 578 CALL wrk_alloc( jpi, jpj, zwx , zwy , zwz ) 579 CALL wrk_alloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) 580 #if defined key_vvl 581 CALL wrk_alloc( jpi, jpj, jpk, ze3f ) 582 #endif 467 CALL wrk_alloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 468 CALL wrk_alloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 583 469 ! 584 470 IF( kt == nit000 ) THEN … … 586 472 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 587 473 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 588 #if ! defined key_vvl589 IF( .NOT.ALLOCATED(ze3f) ) THEN590 ALLOCATE( ze3f(jpi,jpj,jpk) , STAT=ierr )591 IF( lk_mpp ) CALL mpp_sum ( ierr )592 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' )593 ENDIF594 ze3f(:,:,:) = 0._wp595 #endif596 474 ENDIF 597 598 IF( kt == nit000 .OR. lk_vvl ) THEN ! reciprocal of e3 at F-point (masked averaging of e3t over ocean points) 599 600 IF( ln_dynvor_een_old ) THEN ! original formulation 601 DO jk = 1, jpk 602 DO jj = 1, jpjm1 603 DO ji = 1, jpim1 604 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 605 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 606 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = 4.0_wp / ze3 607 END DO 608 END DO 609 END DO 610 ELSE ! new formulation from NEMO 3.6 611 DO jk = 1, jpk 612 DO jj = 1, jpjm1 613 DO ji = 1, jpim1 614 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 615 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 616 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 617 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 618 IF( ze3 /= 0._wp ) ze3f(ji,jj,jk) = zmsk / ze3 619 END DO 620 END DO 621 END DO 622 ENDIF 623 624 CALL lbc_lnk( ze3f, 'F', 1. ) 625 ENDIF 626 627 zfac12 = 1._wp / 12._wp ! Local constant initialization 628 629 630 !CDIR PARALLEL DO PRIVATE( zwx, zwy, zwz, ztnw, ztne, ztsw, ztse ) 475 ! 631 476 ! ! =============== 632 477 DO jk = 1, jpkm1 ! Horizontal slab 633 478 ! ! =============== 634 635 ! Potential vorticity and horizontal fluxes 636 ! ----------------------------------------- 637 SELECT CASE( kvor ) ! vorticity considered 638 CASE ( 1 ) ! planetary vorticity (Coriolis) 639 zwz(:,:) = ff(:,:) * ze3f(:,:,jk) 640 CASE ( 2 ) ! relative vorticity 641 zwz(:,:) = rotn(:,:,jk) * ze3f(:,:,jk) 642 CASE ( 3 ) ! metric term 479 ! 480 SELECT CASE( nn_een_e3f ) ! == reciprocal of e3 at F-point 481 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 482 DO jj = 1, jpjm1 483 DO ji = 1, fs_jpim1 ! vector opt. 484 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 485 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 486 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4.0_wp / ze3 487 ELSE ; z1_e3f(ji,jj) = 0.0_wp 488 ENDIF 489 END DO 490 END DO 491 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 492 DO jj = 1, jpjm1 493 DO ji = 1, fs_jpim1 ! vector opt. 494 ze3 = ( fse3t(ji,jj+1,jk)*tmask(ji,jj+1,jk) + fse3t(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 495 & + fse3t(ji,jj ,jk)*tmask(ji,jj ,jk) + fse3t(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 496 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 497 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 498 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3 499 ELSE ; z1_e3f(ji,jj) = 0.0_wp 500 ENDIF 501 END DO 502 END DO 503 END SELECT 504 ! 505 SELECT CASE( kvor ) !== vorticity considered ==! 506 CASE ( np_COR ) !* Coriolis (planetary vorticity) 507 DO jj = 1, jpjm1 508 DO ji = 1, fs_jpim1 ! vector opt. 509 zwz(ji,jj) = ff(ji,jj) * z1_e3f(ji,jj) 510 END DO 511 END DO 512 CASE ( np_RVO ) !* relative vorticity 513 DO jj = 1, jpjm1 514 DO ji = 1, fs_jpim1 ! vector opt. 515 zwz(ji,jj) = ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 516 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 517 & * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 518 END DO 519 END DO 520 CASE ( np_MET ) !* metric term 643 521 DO jj = 1, jpjm1 644 522 DO ji = 1, fs_jpim1 ! vector opt. 645 523 zwz(ji,jj) = ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 646 524 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 647 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) * ze3f(ji,jj,jk) 648 END DO 649 END DO 650 CALL lbc_lnk( zwz, 'F', 1. ) 651 CASE ( 4 ) ! total (relative + planetary vorticity) 652 zwz(:,:) = ( rotn(:,:,jk) + ff(:,:) ) * ze3f(:,:,jk) 653 CASE ( 5 ) ! total (coriolis + metric) 654 DO jj = 1, jpjm1 655 DO ji = 1, fs_jpim1 ! vector opt. 656 zwz(ji,jj) = ( ff (ji,jj) & 657 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 658 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 659 & * 0.5 / ( e1f(ji,jj) * e2f(ji,jj) ) & 660 & ) * ze3f(ji,jj,jk) 661 END DO 662 END DO 663 CALL lbc_lnk( zwz, 'F', 1. ) 525 & * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 526 END DO 527 END DO 528 CASE ( np_CRV ) !* Coriolis + relative vorticity 529 DO jj = 1, jpjm1 530 DO ji = 1, fs_jpim1 ! vector opt. 531 zwz(ji,jj) = ( ff(ji,jj) + ( e2v(ji+1,jj ) * vn(ji+1,jj ,jk) - e2v(ji,jj) * vn(ji,jj,jk) & 532 & - e1u(ji ,jj+1) * un(ji ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk) ) & 533 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 534 END DO 535 END DO 536 CASE ( np_CME ) !* Coriolis + metric 537 DO jj = 1, jpjm1 538 DO ji = 1, fs_jpim1 ! vector opt. 539 zwz(ji,jj) = ( ff(ji,jj) & 540 & + ( ( vn(ji+1,jj ,jk) + vn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 541 & - ( un(ji ,jj+1,jk) + un (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 542 & * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 543 END DO 544 END DO 545 CASE DEFAULT ! error 546 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 664 547 END SELECT 665 548 ! 549 CALL lbc_lnk( zwz, 'F', 1. ) 550 ! 551 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 552 DO jj = 1, jpjm1 553 DO ji = 1, fs_jpim1 ! vector opt. 554 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 555 END DO 556 END DO 557 ENDIF 558 ! 559 ! !== horizontal fluxes ==! 666 560 zwx(:,:) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 667 561 zwy(:,:) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 668 562 669 ! Compute and add the vorticity term trend 670 ! ---------------------------------------- 563 ! !== compute and add the vorticity term trend =! 671 564 jj = 2 672 565 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 673 DO ji = 2, jpi 566 DO ji = 2, jpi ! split in 2 parts due to vector opt. 674 567 ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 675 568 ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 687 580 DO jj = 2, jpjm1 688 581 DO ji = fs_2, fs_jpim1 ! vector opt. 689 zua = + zfac12 /e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &690 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )691 zva = - zfac12 /e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &692 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )582 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 583 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 584 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 585 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 693 586 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 694 587 pva(ji,jj,jk) = pva(ji,jj,jk) + zva … … 698 591 END DO ! End of slab 699 592 ! ! =============== 700 CALL wrk_dealloc( jpi, jpj, zwx , zwy , zwz ) 701 CALL wrk_dealloc( jpi, jpj, ztnw, ztne, ztsw, ztse ) 702 #if defined key_vvl 703 CALL wrk_dealloc( jpi, jpj, jpk, ze3f ) 704 #endif 593 ! 594 CALL wrk_dealloc( jpi,jpj, zwx , zwy , zwz , z1_e3f ) 595 CALL wrk_dealloc( jpi,jpj, ztnw, ztne, ztsw, ztse ) 705 596 ! 706 597 IF( nn_timing == 1 ) CALL timing_stop('vor_een') … … 720 611 INTEGER :: ios ! Local integer output status for namelist read 721 612 !! 722 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, ln_dynvor_een_old613 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, ln_dynvor_een, nn_een_e3f, ln_dynvor_msk 723 614 !!---------------------------------------------------------------------- 724 615 … … 737 628 WRITE(numout,*) '~~~~~~~~~~~~' 738 629 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 739 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 740 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 741 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 742 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 743 WRITE(numout,*) ' enstrophy and energy conserving scheme (old) ln_dynvor_een_old= ', ln_dynvor_een_old 630 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 631 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 632 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 633 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 634 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 635 WRITE(numout,*) ' masked (=1) or unmasked(=0) vorticity ln_dynvor_msk = ', ln_dynvor_msk 744 636 ENDIF 745 637 638 !!gm this should be removed when choosing a unique strategy for fmask at the coast 746 639 ! If energy, enstrophy or mixed advection of momentum in vector form change the value for masks 747 640 ! at angles with three ocean points and one land point 641 IF(lwp) WRITE(numout,*) 642 IF(lwp) WRITE(numout,*) ' namlbc: change fmask value in the angles (T) ln_vorlat = ', ln_vorlat 748 643 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 749 644 DO jk = 1, jpk … … 759 654 ! 760 655 ENDIF 761 762 ioptio = 0 ! Control of vorticity scheme options 763 IF( ln_dynvor_ene ) ioptio = ioptio + 1 764 IF( ln_dynvor_ens ) ioptio = ioptio + 1 765 IF( ln_dynvor_mix ) ioptio = ioptio + 1 766 IF( ln_dynvor_een ) ioptio = ioptio + 1 767 IF( ln_dynvor_een_old ) ioptio = ioptio + 1 768 IF( lk_esopa ) ioptio = 1 769 656 !!gm end 657 658 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 659 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 660 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 661 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 662 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 663 ! 770 664 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 771 772 ! ! Set nvor (type of scheme for vorticity) 773 IF( ln_dynvor_ene ) nvor = 0 774 IF( ln_dynvor_ens ) nvor = 1 775 IF( ln_dynvor_mix ) nvor = 2 776 IF( ln_dynvor_een .or. ln_dynvor_een_old ) nvor = 3 777 IF( lk_esopa ) nvor = -1 778 779 ! ! Set ncor, nrvm, ntot (type of vorticity) 780 IF(lwp) WRITE(numout,*) 781 ncor = 1 665 ! 666 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) 667 ncor = np_COR 782 668 IF( ln_dynadv_vec ) THEN 783 669 IF(lwp) WRITE(numout,*) ' Vector form advection : vorticity = Coriolis + relative vorticity' 784 nrvm = 2785 ntot = 4670 nrvm = np_RVO ! relative vorticity 671 ntot = np_CRV ! relative + planetary vorticity 786 672 ELSE 787 673 IF(lwp) WRITE(numout,*) ' Flux form advection : vorticity = Coriolis + metric term' 788 nrvm = 3789 ntot = 5674 nrvm = np_MET ! metric term 675 ntot = np_CME ! Coriolis + metric term 790 676 ENDIF 791 677 792 678 IF(lwp) THEN ! Print the choice 793 679 WRITE(numout,*) 794 IF( nvor == 0 ) WRITE(numout,*) ' vorticity scheme : energy conserving scheme' 795 IF( nvor == 1 ) WRITE(numout,*) ' vorticity scheme : enstrophy conserving scheme' 796 IF( nvor == 2 ) WRITE(numout,*) ' vorticity scheme : mixed enstrophy/energy conserving scheme' 797 IF( nvor == 3 ) WRITE(numout,*) ' vorticity scheme : energy and enstrophy conserving scheme' 798 IF( nvor == -1 ) WRITE(numout,*) ' esopa test: use all lateral physics options' 680 IF( nvor_scheme == np_ENE ) WRITE(numout,*) ' vorticity scheme ==>> energy conserving scheme' 681 IF( nvor_scheme == np_ENS ) WRITE(numout,*) ' vorticity scheme ==>> enstrophy conserving scheme' 682 IF( nvor_scheme == np_MIX ) WRITE(numout,*) ' vorticity scheme ==>> mixed enstrophy/energy conserving scheme' 683 IF( nvor_scheme == np_EEN ) WRITE(numout,*) ' vorticity scheme ==>> energy and enstrophy conserving scheme' 799 684 ENDIF 800 685 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r5120 r5836 49 49 !! 50 50 !! ** Method : The now vertical advection of momentum is given by: 51 !! w dz(u) = ua + 1/(e1 u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ]52 !! w dz(v) = va + 1/(e1 v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ]51 !! w dz(u) = ua + 1/(e1e2u*e3u) mk+1[ mi(e1e2t*wn) dk(un) ] 52 !! w dz(v) = va + 1/(e1e2v*e3v) mk+1[ mj(e1e2t*wn) dk(vn) ] 53 53 !! Add this trend to the general trend (ua,va): 54 54 !! (ua,va) = (ua,va) + w dz(u,v) … … 85 85 DO jj = 2, jpj ! vertical fluxes 86 86 DO ji = fs_2, jpi ! vector opt. 87 zww(ji,jj) = 0.25_wp * e1 t(ji,jj) *e2t(ji,jj) * wn(ji,jj,jk)87 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 88 88 END DO 89 89 END DO … … 121 121 DO ji = fs_2, fs_jpim1 ! vector opt. 122 122 ! ! vertical momentum advective trends 123 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )124 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )123 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 124 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 125 125 ! ! add the trends to the general momentum trends 126 126 ua(ji,jj,jk) = ua(ji,jj,jk) + zua … … 146 146 ! 147 147 END SUBROUTINE dyn_zad 148 148 149 149 150 SUBROUTINE dyn_zad_zts ( kt ) … … 182 183 IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts') 183 184 ! 184 CALL wrk_alloc( jpi,jpj,jpk, zwuw, zwvw, zww )185 CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs )185 CALL wrk_alloc( jpi,jpj,jpk, zwuw, zwvw, zww ) 186 CALL wrk_alloc( jpi,jpj,jpk,3, zus , zvs ) 186 187 ! 187 188 IF( kt == nit000 ) THEN … … 205 206 DO jj = 2, jpj 206 207 DO ji = fs_2, jpi ! vector opt. 207 zww(ji,jj,jk) = 0.25_wp * e1 t(ji,jj) *e2t(ji,jj) * wn(ji,jj,jk)208 zww(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 208 209 END DO 209 210 END DO 210 211 END DO 211 ! 212 ! Surface and bottom advective fluxes set to zero 213 DO jj = 2, jpjm1 212 213 DO jj = 2, jpjm1 ! Surface and bottom advective fluxes set to zero 214 214 DO ji = fs_2, fs_jpim1 ! vector opt. 215 !!gm missing ISF boundary condition 215 216 zwuw(ji,jj, 1 ) = 0._wp 216 217 zwvw(ji,jj, 1 ) = 0._wp … … 251 252 DO ji = fs_2, fs_jpim1 ! vector opt. 252 253 ! ! vertical momentum advective trends 253 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1 u(ji,jj) *e2u(ji,jj) * fse3u(ji,jj,jk) )254 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1 v(ji,jj) *e2v(ji,jj) * fse3v(ji,jj,jk) )254 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) ) 255 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1e2v(ji,jj) * fse3v(ji,jj,jk) ) 255 256 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 256 257 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts … … 283 284 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 284 285 ! 285 CALL wrk_dealloc( jpi,jpj,jpk, zwuw, zwvw, zww )286 CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs )286 CALL wrk_dealloc( jpi,jpj,jpk, zwuw, zwvw, zww ) 287 CALL wrk_dealloc( jpi,jpj,jpk,3, zus , zvs ) 287 288 ! 288 289 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts') -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r4990 r5836 9 9 10 10 !!---------------------------------------------------------------------- 11 !! dyn_zdf : Update the momentum trend with the vertical diffusion12 !! dyn_zdf_init : initializations of the vertical diffusion scheme11 !! dyn_zdf : Update the momentum trend with the vertical diffusion 12 !! dyn_zdf_init : initializations of the vertical diffusion scheme 13 13 !!---------------------------------------------------------------------- 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physics variables 17 18 USE dynzdf_exp ! vertical diffusion: explicit (dyn_zdf_exp routine) 19 USE dynzdf_imp ! vertical diffusion: implicit (dyn_zdf_imp routine) 20 21 USE ldfdyn_oce ! ocean dynamics: lateral physics 22 USE trd_oce ! trends: ocean variables 23 USE trddyn ! trend manager: dynamics 24 USE in_out_manager ! I/O manager 25 USE lib_mpp ! MPP library 26 USE prtctl ! Print control 27 USE wrk_nemo ! Memory Allocation 28 USE timing ! Timing 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physics variables 17 USE dynzdf_exp ! vertical diffusion: explicit (dyn_zdf_exp routine) 18 USE dynzdf_imp ! vertical diffusion: implicit (dyn_zdf_imp routine) 19 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 20 USE trd_oce ! trends: ocean variables 21 USE trddyn ! trend manager: dynamics 22 ! 23 USE in_out_manager ! I/O manager 24 USE lib_mpp ! MPP library 25 USE prtctl ! Print control 26 USE wrk_nemo ! Memory Allocation 27 USE timing ! Timing 29 28 30 29 IMPLICIT NONE … … 61 60 !!--------------------------------------------------------------------- 62 61 ! 63 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf')62 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf') 64 63 ! 65 64 ! ! set time step … … 79 78 CASE ( 1 ) ; CALL dyn_zdf_imp( kt, r2dt ) ! implicit scheme 80 79 ! 81 CASE ( -1 ) ! esopa: test all possibility with control print82 CALL dyn_zdf_exp( kt, r2dt )83 CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf0 - Ua: ', mask1=umask, &84 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )85 CALL dyn_zdf_imp( kt, r2dt )86 CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf1 - Ua: ', mask1=umask, &87 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )88 80 END SELECT 89 81 … … 96 88 ! ! print mean trends (used for debugging) 97 89 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zdf - Ua: ', mask1=umask, & 98 &tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )99 !100 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf')90 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 91 ! 92 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf') 101 93 ! 102 94 END SUBROUTINE dyn_zdf … … 114 106 USE zdftke 115 107 USE zdfgls 116 USE zdfkpp117 108 !!---------------------------------------------------------------------- 118 109 ! … … 123 114 ! 124 115 ! Force implicit schemes 125 IF( lk_zdftke .OR. lk_zdfgls .OR. lk_zdfkpp ) nzdf = 1 ! TKE, GLS or KPP physics 126 IF( ln_dynldf_iso ) nzdf = 1 ! iso-neutral lateral physics 127 IF( ln_dynldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate 128 ! 129 IF( lk_esopa ) nzdf = -1 ! Esopa key: All schemes used 116 IF( lk_zdftke .OR. lk_zdfgls ) nzdf = 1 ! TKE or GLS physics 117 IF( ln_dynldf_iso ) nzdf = 1 ! iso-neutral lateral physics 118 IF( ln_dynldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate 130 119 ! 131 120 IF(lwp) THEN ! Print the choice … … 133 122 WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme' 134 123 WRITE(numout,*) '~~~~~~~~~~~' 135 IF( nzdf == -1 ) WRITE(numout,*) ' ESOPA test All scheme used'136 124 IF( nzdf == 0 ) WRITE(numout,*) ' Explicit time-splitting scheme' 137 125 IF( nzdf == 1 ) WRITE(numout,*) ' Implicit (euler backward) scheme' -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r5120 r5836 209 209 !----------------------------------------------------------------------- 210 210 ! 211 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 212 DO jk = 2, jpkm1 211 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 213 212 DO jj = 2, jpjm1 214 213 DO ji = fs_2, fs_jpim1 ! vector opt. … … 309 308 !----------------------------------------------------------------------- 310 309 ! 311 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 312 DO jk = 2, jpkm1 310 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 313 311 DO jj = 2, jpjm1 314 312 DO ji = fs_2, fs_jpim1 ! vector opt. -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5656 r5836 20 20 USE sbc_oce ! surface boundary condition: ocean 21 21 USE domvvl ! Variable volume 22 USE divcur ! hor. divergence and curl (div & cur routines) 23 USE restart ! only for lrst_oce 24 USE in_out_manager ! I/O manager 25 USE prtctl ! Print control 26 USE phycst 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 28 USE lib_mpp ! MPP library 22 USE divhor ! horizontal divergence 23 USE phycst ! physical constants 29 24 USE bdy_oce 30 25 USE bdy_par … … 36 31 USE asminc ! Assimilation increment 37 32 #endif 33 USE in_out_manager ! I/O manager 34 USE restart ! only for lrst_oce 35 USE prtctl ! Print control 36 USE lbclnk ! ocean lateral boundary condition (or mpp link) 37 USE lib_mpp ! MPP library 38 38 USE wrk_nemo ! Memory Allocation 39 39 USE timing ! Timing … … 66 66 !! by the time step. 67 67 !! 68 !! ** action : ssha :after sea surface height68 !! ** action : ssha, after sea surface height 69 69 !! 70 70 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 71 71 !!---------------------------------------------------------------------- 72 ! 73 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 74 INTEGER, INTENT(in) :: kt ! time step 72 INTEGER, INTENT(in) :: kt ! time step 75 73 ! 76 INTEGER :: jk ! dummy loop indice 77 REAL(wp) :: z2dt, z1_rau0 ! local scalars 78 !!---------------------------------------------------------------------- 79 ! 80 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 81 ! 82 CALL wrk_alloc( jpi, jpj, zhdiv ) 74 INTEGER :: jk ! dummy loop indice 75 REAL(wp) :: z2dt, zcoef ! local scalars 76 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv ! 2D workspace 77 !!---------------------------------------------------------------------- 78 ! 79 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 80 ! 81 CALL wrk_alloc( jpi,jpj, zhdiv ) 83 82 ! 84 83 IF( kt == nit000 ) THEN 85 !86 84 IF(lwp) WRITE(numout,*) 87 85 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 88 86 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 89 ! 90 ENDIF 91 ! 92 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 87 ENDIF 88 ! 89 CALL div_hor( kt ) ! Horizontal divergence 93 90 ! 94 91 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) … … 106 103 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 107 104 ! 108 z 1_rau0= 0.5_wp * r1_rau0109 ssha(:,:) = ( sshb(:,:) - z2dt * ( z 1_rau0* ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)105 zcoef = 0.5_wp * r1_rau0 106 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 110 107 111 108 #if ! defined key_dynspg_ts 112 109 ! These lines are not necessary with time splitting since 113 110 ! boundary condition on sea level is set during ts loop 114 # if defined key_agrif111 # if defined key_agrif 115 112 CALL agrif_ssh( kt ) 116 # endif117 # if defined key_bdy118 IF (lk_bdy) THEN119 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary120 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries121 ENDIF 122 # endif113 # endif 114 # if defined key_bdy 115 IF( lk_bdy ) THEN 116 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 117 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 118 ENDIF 119 # endif 123 120 #endif 124 121 125 122 #if defined key_asminc 126 ! ! Include the IAU weighted SSH increment 127 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 123 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN ! Include the IAU weighted SSH increment 128 124 CALL ssh_asm_inc( kt ) 129 125 ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 130 126 ENDIF 131 127 #endif 132 133 128 ! !------------------------------! 134 129 ! ! outputs ! … … 159 154 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 160 155 !!---------------------------------------------------------------------- 161 ! 162 INTEGER, INTENT(in) :: kt ! time step 156 INTEGER, INTENT(in) :: kt ! time step 157 ! 158 INTEGER :: ji, jj, jk ! dummy loop indices 159 REAL(wp) :: z1_2dt ! local scalars 163 160 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 164 161 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 165 ! 166 INTEGER :: ji, jj, jk ! dummy loop indices 167 REAL(wp) :: z1_2dt ! local scalars 168 !!---------------------------------------------------------------------- 169 170 IF( nn_timing == 1 ) CALL timing_start('wzv') 162 !!---------------------------------------------------------------------- 163 ! 164 IF( nn_timing == 1 ) CALL timing_start('wzv') 171 165 ! 172 166 IF( kt == nit000 ) THEN 173 !174 167 IF(lwp) WRITE(numout,*) 175 168 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' … … 177 170 ! 178 171 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 179 !180 172 ENDIF 181 173 ! !------------------------------! … … 193 185 DO jj = 2, jpjm1 194 186 DO ji = fs_2, fs_jpim1 ! vector opt. 195 zhdiv(ji,jj,jk) = r1_e1 2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )187 zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 196 188 END DO 197 189 END DO … … 216 208 217 209 #if defined key_bdy 218 IF (lk_bdy) THEN210 IF( lk_bdy ) THEN 219 211 DO jk = 1, jpkm1 220 212 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) … … 224 216 ! 225 217 IF( nn_timing == 1 ) CALL timing_stop('wzv') 226 227 218 ! 228 219 END SUBROUTINE wzv 220 229 221 230 222 SUBROUTINE ssh_swp( kt ) … … 265 257 sshb(:,:) = sshn(:,:) ! before <-- now 266 258 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 259 ! 267 260 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 268 261 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered
Note: See TracChangeset
for help on using the changeset viewer.