Changeset 13154
- Timestamp:
- 2020-06-24T15:31:32+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE
- Files:
-
- 4 edited
-
dynadv.F90 (modified) (1 diff)
-
dynkeg.F90 (modified) (2 diffs)
-
dynvor.F90 (modified) (15 diffs)
-
step.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynadv.F90
r13151 r13154 74 74 CASE( np_VEC_c2 ) 75 75 CALL dyn_keg ( kt, nn_dynkeg, Kmm, puu, pvv, Krhs ) ! vector form : horizontal gradient of kinetic energy 76 !!an SWE : w = 0 76 77 CALL dyn_zad ( kt, Kmm, puu, pvv, Krhs ) ! vector form : vertical advection 77 78 CASE( np_FLX_c2 ) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
r13151 r13154 29 29 30 30 PUBLIC dyn_keg ! routine called by step module 31 PUBLIC dyn_kegAD ! routine called by step module 31 32 32 33 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) … … 155 156 ! 156 157 END SUBROUTINE dyn_keg 157 158 159 160 SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs ) 161 !!---------------------------------------------------------------------- 162 !! *** ROUTINE dyn_kegAD *** 163 !! 164 !! ** Purpose : Compute the now momentum trend due to the horizontal 165 !! gradient of the horizontal kinetic energy and add it to the 166 !! general momentum trend. 167 !! 168 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 169 !! conserve kinetic energy. Compute the now horizontal kinetic energy 170 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 171 !! * kscheme = nkeg_HW : Hollingsworth correction following 172 !! Arakawa (2001). The now horizontal kinetic energy is given by: 173 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 ) 174 !! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ] 175 !! 176 !! Take its horizontal gradient and add it to the general momentum 177 !! trend. 178 !! u(rhs) = u(rhs) - 1/e1u di[ zhke ] 179 !! v(rhs) = v(rhs) - 1/e2v dj[ zhke ] 180 !! 181 !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend 182 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 183 !! 184 !! ** References : Arakawa, A., International Geophysics 2001. 185 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 186 !!---------------------------------------------------------------------- 187 ! 188 INTEGER , INTENT( in ) :: kt ! ocean time-step index 189 INTEGER , INTENT( in ) :: kscheme ! =0/1/2 type of KEG scheme 190 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: puu, pvv ! ocean velocities at Kmm 191 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! RHS 192 ! 193 INTEGER :: ji, jj, jk ! dummy loop indices 194 REAL(wp) :: zu, zv ! local scalars 195 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke 196 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 197 !!---------------------------------------------------------------------- 198 ! 199 IF( ln_timing ) CALL timing_start('dyn_keg') 200 ! 201 IF( kt == nit000 ) THEN 202 IF(lwp) WRITE(numout,*) 203 IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme 204 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 205 ENDIF 206 207 zhke(:,:,jpk) = 0._wp 208 209 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 210 !!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2) 211 CASE ( nkeg_C2_wpo ) !-- Standard scheme --! 212 DO_3D_01_01( 1, jpkm1 ) 213 zu = ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 214 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) 215 zv = ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 216 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) 217 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 218 END_3D 219 !!an45 220 ! 221 CASE ( nkeg_C2 ) !-- Standard scheme --! 222 DO_3D_01_01( 1, jpkm1 ) 223 zu = puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 224 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) 225 zv = pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 226 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) 227 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 228 END_3D 229 !!an 00_00 ? 230 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 231 DO_3D_00_00( 1, jpkm1 ) 232 zu = 8._wp * ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 233 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) & 234 & + ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) & 235 & + ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) ) * ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) ) 236 ! 237 zv = 8._wp * ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 238 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) & 239 & + ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) & 240 & + ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) ) * ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) ) 241 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 242 END_3D 243 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 244 ! 245 END SELECT 246 ! 247 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 248 ! 249 DO_3D_00_00( 1, jpkm1 ) 250 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 251 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 252 END_3D 253 ! 254 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 255 ! 256 DO_3D_00_00( 1, jpkm1 ) 257 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 258 END_3D 259 ! 260 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 261 ! 262 DO_3D_00_00( 1, jpkm1 ) 263 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 264 END_3D 265 ! 266 ENDIF 267 IF( ln_timing ) CALL timing_stop('dyn_kegAD') 268 ! 269 END SUBROUTINE dyn_kegAD 158 270 !!====================================================================== 159 271 END MODULE dynkeg -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynvor.F90
r13151 r13154 22 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 23 23 !! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity 24 !! 4.x ! 2020-03 (G. Madec, A. Nasser) alternate direction computation of vorticity tendancy 25 !! ! for ENS, ENE 24 26 !!---------------------------------------------------------------------- 25 27 … … 45 47 USE lib_mpp ! MPP library 46 48 USE timing ! Timing 49 !!anAD only 50 USE dynkeg, ONLY : dyn_kegAD 51 USE dynadv, ONLY : nn_dynkeg 47 52 48 53 IMPLICIT NONE … … 53 58 54 59 ! !!* Namelist namdyn_vor: vorticity term 55 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 56 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 57 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 58 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 59 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 60 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 61 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 62 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 60 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 61 LOGICAL, PUBLIC :: ln_dynvor_ens_adVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 62 LOGICAL, PUBLIC :: ln_dynvor_ens_adKE = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 63 LOGICAL, PUBLIC :: ln_dynvor_ens_adKEVO = .FALSE. !: AD enstrophy conserving scheme (ENS_ad) 64 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 65 LOGICAL, PUBLIC :: ln_dynvor_ene_adVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 66 LOGICAL, PUBLIC :: ln_dynvor_ene_adKE = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 67 LOGICAL, PUBLIC :: ln_dynvor_ene_adKEVO = .FALSE. !: f-point AD energy conserving scheme (ENE_ad) 68 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 69 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 70 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 71 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 72 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 73 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 63 74 64 75 INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme … … 70 81 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 71 82 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 72 73 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 83 !!an 84 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKE = 11 ! ENS scheme - AD scheme (KE only) 85 INTEGER, PUBLIC, PARAMETER :: np_ENS_adVO = 12 ! ENS scheme - AD scheme (VOR only) 86 INTEGER, PUBLIC, PARAMETER :: np_ENS_adKEVO = 13 ! ENS scheme - AD scheme (KE+VOR) 87 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKE = 21 ! ENE scheme - AD scheme (KE only) 88 INTEGER, PUBLIC, PARAMETER :: np_ENE_adVO = 22 ! ENE scheme - AD scheme (VOR only) 89 INTEGER, PUBLIC, PARAMETER :: np_ENE_adKEVO = 23 ! ENE scheme - AD scheme (KE+VOR) 90 !!an 91 92 !!an ds step on pourra spécifier la valeur de ntot = np_COR ou np_COR + np_RVO 93 INTEGER, PUBLIC :: ncor, nrvm, ntot ! choice of calculated vorticity 74 94 ! ! associated indices: 75 95 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) … … 99 119 CONTAINS 100 120 101 SUBROUTINE dyn_vor( kt, K mm, puu, pvv, Krhs)121 SUBROUTINE dyn_vor( kt, Kbb, Kmm, puu, pvv, Krhs) 102 122 !!---------------------------------------------------------------------- 103 123 !! … … 109 129 !! for futher diagnostics (l_trddyn=T) 110 130 !!---------------------------------------------------------------------- 111 INTEGER , INTENT( in ) :: kt ! ocean time-step index 112 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 113 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 131 INTEGER :: ji, jj, jk ! dummy loop indice 132 INTEGER , INTENT( in ) :: kt ! ocean time-step index 133 INTEGER , INTENT( in ) :: Kmm, Krhs, Kbb ! ocean time level indices 134 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 114 135 ! 115 136 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 137 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zuu, zvv 116 138 !!---------------------------------------------------------------------- 117 139 ! … … 167 189 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 168 190 CASE( np_ENE ) !* energy conserving scheme 191 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 192 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 193 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 194 195 CASE( np_ENE_adVO ) !* energy conserving scheme with alternating direction 196 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 197 198 !== Alternative direction - VOR only ==! 199 200 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 201 202 ALLOCATE( zuu(jpi,jpj,jpk) ) 203 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 204 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 205 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 206 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 207 DEALLOCATE( zuu ) 208 ELSE 209 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 210 211 ALLOCATE( zvv(jpi,jpj,jpk) ) 212 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 213 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 214 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 215 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 216 DEALLOCATE( zvv ) 217 ENDIF 218 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 219 CASE( np_ENE_adKE ) !* energy conserving scheme with alternating direction 220 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 221 222 !== Alternative direction - KEG only ==! 223 169 224 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 225 226 ALLOCATE( zuu(jpi,jpj,jpk) ) 227 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 228 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 229 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 230 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 231 DEALLOCATE( zuu ) 232 ELSE 233 234 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 235 236 ALLOCATE( zvv(jpi,jpj,jpk) ) 237 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 238 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 239 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 240 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 241 DEALLOCATE( zvv ) 242 ENDIF 170 243 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 244 245 CASE( np_ENE_adKEVO ) !* energy conserving scheme with alternating direction 246 ! equivalent ADE pas sur coriolis 247 ! appel vor_ene(ncor) 248 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 249 250 !== Alternative direction - KE + VOR ==! 251 252 ALLOCATE( zuu(jpi,jpj,jpk) ) 253 !puis call que sur nrvm 254 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 255 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 256 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 257 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 258 !puis call que sur nrvm 259 CALL vor_ene( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 260 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 261 DEALLOCATE( zuu ) 262 ELSE 263 264 ALLOCATE( zvv(jpi,jpj,jpk) ) 265 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 266 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 267 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 268 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 269 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 270 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 271 DEALLOCATE( zvv ) 272 ENDIF 171 273 CASE( np_ENS ) !* enstrophy conserving scheme 274 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 275 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 276 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 277 CASE( np_ENS_adVO ) !* enstrophy conserving scheme with alternating direction 278 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 279 280 !== Alternative direction - VOR only ==! 281 282 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 283 284 ALLOCATE( zuu(jpi,jpj,jpk) ) 285 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 286 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 287 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 288 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 289 DEALLOCATE( zuu ) 290 ELSE 291 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs), pv_rhs=pvv(:,:,:,Krhs) ) 292 293 ALLOCATE( zvv(jpi,jpj,jpk) ) 294 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 295 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 296 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 297 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , zvv(:,:,:), pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 298 DEALLOCATE( zvv ) 299 ENDIF 300 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 301 CASE( np_ENS_adKE ) !* enstrophy conserving scheme with alternating direction 302 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 303 304 !== Alternative direction - KEG only ==! 305 172 306 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 307 308 ALLOCATE( zuu(jpi,jpj,jpk) ) 309 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 310 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 311 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 312 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 313 DEALLOCATE( zuu ) 314 ELSE 315 316 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 317 318 ALLOCATE( zvv(jpi,jpj,jpk) ) 319 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 320 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 321 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 322 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 323 DEALLOCATE( zvv ) 324 ENDIF 325 CASE( np_ENS_adKEVO ) !* enstrophy conserving scheme with alternating direction 326 IF( MOD( kt , 2 ) == 1 ) THEN ! even time step: u-vor then v-vor components 327 328 !== Alternative direction - KE + VOR ==! 329 330 ALLOCATE( zuu(jpi,jpj,jpk) ) 331 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 332 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pu_rhs=puu(:,:,:,Krhs) ) ! 333 zuu(:,:,:) = puu(:,:,:,Kbb) + rDt * puu(:,:,:,Krhs) * umask(:,:,:) 334 CALL lbc_lnk( 'dynvor', zuu(:,:,:) , 'U', -1._wp ) 335 CALL vor_ens( kt, Kmm, ntot, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 336 CALL dyn_kegAD( kt, nn_dynkeg, zuu(:,:,:) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 337 DEALLOCATE( zuu ) 338 ELSE 339 340 ALLOCATE( zvv(jpi,jpj,jpk) ) 341 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) ! compute and add vv-vorticity trend 342 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , pv_rhs=pvv(:,:,:,Krhs) ) 343 zvv(:,:,:) = pvv(:,:,:,Kbb) + rDt * pvv(:,:,:,Krhs) * vmask(:,:,:) 344 CALL lbc_lnk( 'dynvor', zvv(:,:,:) , 'V', -1._wp ) 345 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) ! compute and add uu-vorticity trend 346 CALL dyn_kegAD( kt, nn_dynkeg, puu(:,:,:,Kmm), zvv(:,:,:) , pu_rhs=puu(:,:,:,Krhs) ) 347 DEALLOCATE( zvv ) 348 ENDIF 173 349 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 174 350 CASE( np_MIX ) !* mixed ene-ens scheme … … 319 495 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 320 496 !!---------------------------------------------------------------------- 321 INTEGER , INTENT(in ) :: kt ! ocean time-step index322 INTEGER , INTENT(in ) :: Kmm ! ocean time level index323 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric324 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities325 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend497 INTEGER , INTENT(in ) :: kt ! ocean time-step index 498 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 499 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 500 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 501 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 326 502 ! 327 503 INTEGER :: ji, jj, jk ! dummy loop indices 328 504 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 329 505 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! 2D workspace 506 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t ! unpenalised scale factors 330 507 !!---------------------------------------------------------------------- 331 508 ! … … 377 554 END SELECT 378 555 ! 379 ! !== horizontal fluxes and potential vorticity ==! 380 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 381 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 382 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 383 ! 384 ! !== compute and add the vorticity term trend =! 385 DO_2D_00_00 386 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 387 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 388 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 389 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 390 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 391 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 392 END_2D 556 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 557 ! 558 ! !== horizontal fluxes and potential vorticity ==! 559 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 560 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 561 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 562 ! 563 ! !== compute and add the vorticity term trend =! 564 DO_2D_00_00 565 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 566 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 567 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 568 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 569 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 570 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 571 END_2D 572 ! 573 ! 574 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 575 ! 576 ! 577 ! !== horizontal fluxes and potential vorticity ==! 578 !!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée 579 ! on transportera (en AD) de la bonne manière 580 !zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 581 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * vv(:,:,jk,Kmm) 582 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 583 ! 584 ! !== compute and add the vorticity term trend =! 585 DO_2D_00_00 586 zy1 = zwy(ji,jj-1) + zwy(ji+1,jj-1) 587 zy2 = zwy(ji,jj ) + zwy(ji+1,jj ) 588 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 589 END_2D 590 ! 591 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 592 ! 593 ! 594 ! !== horizontal fluxes and potential vorticity ==! 595 !!anAD remplacer pu par uu(Nnn) et on fait correctement la direction alternée 596 !zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 597 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * uu(:,:,jk,Kmm) 598 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 599 ! 600 ! !== compute and add the vorticity term trend =! 601 DO_2D_00_00 602 zx1 = zwx(ji-1,jj) + zwx(ji-1,jj+1) 603 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 604 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 605 END_2D 606 ! 607 ENDIF 393 608 ! ! =============== 394 609 END DO ! End of slab … … 417 632 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. 418 633 !!---------------------------------------------------------------------- 419 INTEGER , INTENT(in ) :: kt ! ocean time-step index420 INTEGER , INTENT(in ) :: Kmm ! ocean time level index421 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric422 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities423 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend634 INTEGER , INTENT(in ) :: kt ! ocean time-step index 635 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 636 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 637 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: pu, pv ! now velocities 638 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 424 639 ! 425 640 INTEGER :: ji, jj, jk ! dummy loop indices … … 475 690 ! 476 691 ! 477 ! !== horizontal fluxes and potential vorticity ==! 478 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 479 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 480 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 481 ! 482 ! !== compute and add the vorticity term trend =! 483 DO_2D_00_00 484 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 485 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 486 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 487 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 488 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 489 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 490 END_2D 692 !!an wut ? v et u 693 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 694 ! 695 ! !== horizontal fluxes and potential vorticity ==! 696 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 697 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 698 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 699 ! 700 ! !== compute and add the vorticity term trend =! 701 DO_2D_00_00 702 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 703 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 704 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 705 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 706 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 707 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 708 END_2D 709 ! 710 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 711 ! 712 ! 713 ! !== horizontal fluxes and potential vorticity ==! 714 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 715 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 716 ! 717 ! !== compute and add the vorticity term trend =! 718 DO_2D_00_00 719 zuav = r1_8 * r1_e1u(ji,jj) * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 720 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) 721 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zuav * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 722 END_2D 723 ! 724 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 725 ! 726 ! 727 ! !== horizontal fluxes and potential vorticity ==! 728 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 729 zwz(:,:) = zwz(:,:) / e3f(:,:,jk) 730 ! 731 ! !== compute and add the vorticity term trend =! 732 DO_2D_00_00 733 zvau =-r1_8 * r1_e2v(ji,jj) * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 734 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) 735 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zvau * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 736 END_2D 737 ! 738 ENDIF 491 739 ! ! =============== 492 740 END DO ! End of slab … … 772 1020 !! 773 1021 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 774 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 1022 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk, & 1023 & ln_dynvor_ens_adVO, ln_dynvor_ens_adKE, ln_dynvor_ens_adKEVO, & ! Alternative direction parameters 1024 & ln_dynvor_ene_adVO, ln_dynvor_ene_adKE, ln_dynvor_ene_adKEVO 775 1025 !!---------------------------------------------------------------------- 776 1026 ! … … 809 1059 DO_3D_10_10( 1, jpk ) 810 1060 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 811 & + tmask(ji,jj ,jk) + tmask(ji+1,jj +1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp1061 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 812 1062 END_3D 813 1063 ! … … 819 1069 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 820 1070 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 1071 IF( ln_dynvor_ens_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adVO ; ENDIF 1072 IF( ln_dynvor_ens_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKE ; ENDIF 1073 IF( ln_dynvor_ens_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS_adKEVO ; ENDIF 821 1074 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 1075 IF( ln_dynvor_ene_adVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adVO ; ENDIF 1076 IF( ln_dynvor_ene_adKE ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKE ; ENDIF 1077 IF( ln_dynvor_ene_adKEVO ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE_adKEVO ; ENDIF 822 1078 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 823 1079 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF … … 866 1122 WRITE(numout,*) 867 1123 SELECT CASE( nvor_scheme ) 868 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 869 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 870 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 871 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 872 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 873 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 1124 1125 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 1126 CASE( np_ENS_adVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adVO) on vorticity only' 1127 CASE( np_ENS_adKE ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKE) on kinetic energy only' 1128 CASE( np_ENS_adKEVO ) ; WRITE(numout,*) ' ==>>> AD enstrophy conserving scheme (ENS_adKEVO) on kinetic energy and vorticity' 1129 1130 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 1131 CASE( np_ENE_adVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adVO) on vorticity only' 1132 CASE( np_ENE_adKE ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKE) on kinetic energy only' 1133 CASE( np_ENE_adKEVO ) ; WRITE(numout,*) ' ==>>> AD energy conserving scheme (Coriolis at F-points) (ENE_adKEVO) on kinetic energy and vorticity' 1134 1135 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 1136 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 1137 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 1138 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 874 1139 END SELECT 875 1140 ENDIF -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/step.F90
r13151 r13154 17 17 USE phycst ! physical constants 18 18 USE usrdef_nam 19 !19 USE lib_mpp ! MPP library 20 20 USE iom ! xIOs server 21 21 … … 66 66 REAL(wp) :: zue3a, zue3n, zue3b ! local scalars 67 67 REAL(wp) :: zve3a, zve3n, zve3b ! - - 68 REAL(wp) :: ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 68 REAL(wp) :: ze3t_tf, ze3u_tf, ze3v_tf, zua, zva 69 REAL(wp), DIMENSION(jpi,jpj) :: zssh 70 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3u, ze3v 69 71 !! --------------------------------------------------------------------- 70 72 #if defined key_agrif … … 141 143 & CALL Agrif_Sponge_dyn ! momentum sponge 142 144 #endif 143 CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS144 145 CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS146 147 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing148 145 149 146 !!an - calcul du gradient de pression horizontal (explicit) 147 zssh(:,:) = ssh(:,:,Nnn) 148 # if defined key_bvp 149 ! gradient and divergence not penalised 150 zssh(:,:) = zssh(:,:) * r1_rpo(:,:,1) 151 #endif 150 152 DO_3D_00_00( 1, jpkm1 ) 151 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)152 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)153 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( zssh(ji+1,jj) - zssh(ji,jj) ) * r1_e1u(ji,jj) 154 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( zssh(ji,jj+1) - zssh(ji,jj) ) * r1_e2v(ji,jj) 153 155 END_3D 154 156 ! 155 ! add wind stress forcing and layer linear friction to the RHS 157 158 ! IF( kstp == nit000 .AND. lwp ) THEN 159 ! WRITE(numout,*) 160 ! WRITE(numout,*) 'step.F90 : classic script used' 161 ! WRITE(numout,*) '~~~~~~~' 162 ! IF( ln_dynvor_ens_adVO .OR. ln_dynvor_ens_adKE .OR. ln_dynvor_ens_adKEVO & 163 ! & .OR. ln_dynvor_ene_adVO .OR. ln_dynvor_ene_adKE .OR. ln_dynvor_ene_adKEVO ) THEN 164 ! CALL ctl_stop('STOP','step : alternative direction asked but classis step' ) 165 ! ENDIF 166 ! ENDIF 167 !!an 168 ! CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS 169 ! 170 ! CALL dyn_vor( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS 171 ! 172 !!an In dynvor, dynkegAD is called even if not AD, so we keep the same step.F90 173 174 CALL dyn_vor( kstp, Nbb, Nnn , uu, vv, Nrhs) ! vorticity ==> RHS 175 176 CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing 177 178 ! add wind stress forcing and layer linear friction to the RHS 179 ze3u(:,:,:) = e3u(:,:,:,Nnn) 180 ze3v(:,:,:) = e3v(:,:,:,Nnn) 181 # if defined key_bvp 182 !ze3u(:,:,:) = ze3u(:,:,:) * r1_rpo(:,:,:) 183 !ze3v(:,:,:) = ze3v(:,:,:) * r1_rpo(:,:,:) 184 #endif 156 185 z1_2rho0 = 0.5_wp * r1_rho0 157 186 DO_3D_00_00(1,jpkm1) 158 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) &187 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / ze3u(ji,jj,jk) & 159 188 & - rn_rfr * uu(ji,jj,jk,Nbb) 160 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) &161 & - rn_rfr * vv(ji,jj,jk,Nbb) 189 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ze3v(ji,jj,jk) & 190 & - rn_rfr * vv(ji,jj,jk,Nbb) 162 191 END_3D 163 !!an 192 ! 193 # if defined key_bvp 194 ! Add frictionnal term - sigma * u 195 ! can be done via sigmaT (mean_ij) 196 DO_3D_00_00(1,jpkm1) 197 uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - bmpu(ji,jj,jk) * uu(ji,jj,jk,Nbb) 198 vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - bmpv(ji,jj,jk) * vv(ji,jj,jk,Nbb) 199 END_3D 200 #endif 201 ! 164 202 165 203 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Note: See TracChangeset
for help on using the changeset viewer.
