Changeset 9528 for branches/2017/dev_merge_2017/NEMOGCM
- Timestamp:
- 2018-04-30T15:39:18+02:00 (6 years ago)
- Location:
- branches/2017/dev_merge_2017/NEMOGCM
- Files:
-
- 5 edited
-
CONFIG/GYRE_BFM/EXP00/namelist_cfg (modified) (1 diff)
-
CONFIG/GYRE_PISCES/EXP00/namelist_cfg (modified) (1 diff)
-
CONFIG/SHARED/namelist_ref (modified) (4 diffs)
-
NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (28 diffs)
-
NEMO/OPA_SRC/DYN/dynvor.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/CONFIG/GYRE_BFM/EXP00/namelist_cfg
r9527 r9528 167 167 &namdyn_vor ! Vorticity / Coriolis scheme (default: NO selection) 168 168 !----------------------------------------------------------------------- 169 ln_dynvor_ene = .true. ! en strophy conserving scheme169 ln_dynvor_ene = .true. ! energy conserving scheme 170 170 / 171 171 !----------------------------------------------------------------------- -
branches/2017/dev_merge_2017/NEMOGCM/CONFIG/GYRE_PISCES/EXP00/namelist_cfg
r9527 r9528 168 168 &namdyn_vor ! Vorticity / Coriolis scheme (default: NO selection) 169 169 !----------------------------------------------------------------------- 170 ln_dynvor_ene = .true. ! en strophy conserving scheme170 ln_dynvor_ene = .true. ! energy conserving scheme 171 171 / 172 172 !----------------------------------------------------------------------- -
branches/2017/dev_merge_2017/NEMOGCM/CONFIG/SHARED/namelist_ref
r9527 r9528 2 2 !! NEMO/OCE : Reference namelist_ref !! 3 3 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 4 !! NEMO/OPA : 1 - run manager (namrun) 5 !! namelists 2 - Domain (namcfg, namdom, namtsd, namcrs, namc1d, namc1d_uvd) 6 !! 3 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl, 4 !! NEMO/OPA : 1 - Domain & run manager (namrun, namcfg, namdom, namtsd, namcrs, namc1d, namc1d_uvd) 5 !! namelists 2 - Surface boundary (namsbc, namsbc_flx, namsbc_blk, namsbc_cpl, 7 6 !! namsbc_sas, namtra_qsr, namsbc_rnf, 8 7 !! namsbc_isf, namsbc_iscpl, namsbc_apr, 9 8 !! namsbc_ssr, namsbc_wave, namberg) 10 !! 4- lateral boundary (namlbc, namagrif, nambdy, nambdy_tide)11 !! 5- top/bot boundary (namdrg, namdrg_top, namdrg_bot, nambbc, nambbl)12 !! 6- Tracer (nameos, namtra_adv, namtra_ldf, namtra_eiv, namtra_dmp)13 !! 7- dynamics (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf)14 !! 8- Vertical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_iwm)15 !! 9 - miscellaneous (nammpp, namctl)16 !! 10 - diagnostics (namnc4, namtrd, namspr, namflo, namhsb, namsto)17 !! 1 1 - Obs & Assim (namobs, nam_asminc)9 !! 3 - lateral boundary (namlbc, namagrif, nambdy, nambdy_tide) 10 !! 4 - top/bot boundary (namdrg, namdrg_top, namdrg_bot, nambbc, nambbl) 11 !! 5 - Tracer (nameos, namtra_adv, namtra_ldf, namtra_eiv, namtra_dmp) 12 !! 6 - dynamics (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 13 !! 7 - Vertical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_gls, namzdf_iwm) 14 !! 8 - diagnostics (namnc4, namtrd, namspr, namflo, namhsb) 15 !! 9 - Obs & Assim (namobs, nam_asminc) 16 !! 10 - miscellaneous (nammpp, namctl, namsto) 18 17 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 19 18 … … 651 650 &namdrg ! top/bottom drag coefficient (default: NO selection) 652 651 !----------------------------------------------------------------------- 653 ln_OFF = .false.! free-slip : Cd = 0 (F => fill namdrg_bot654 ln_lin = .false.! linear drag: Cd = Cd0 Uc0 & namdrg_top)655 ln_non_lin = .false.! non-linear drag: Cd = Cd0 |U|656 ln_loglayer = .false.! logarithmic drag: Cd = vkarmn/log(z/z0) |U|657 ! 658 ln_drgimp = .true.! implicit top/bottom friction flag652 ln_OFF = .false. ! free-slip : Cd = 0 (F => fill namdrg_bot 653 ln_lin = .false. ! linear drag: Cd = Cd0 Uc0 & namdrg_top) 654 ln_non_lin = .false. ! non-linear drag: Cd = Cd0 |U| 655 ln_loglayer = .false. ! logarithmic drag: Cd = vkarmn/log(z/z0) |U| 656 ! 657 ln_drgimp = .true. ! implicit top/bottom friction flag 659 658 / 660 659 !----------------------------------------------------------------------- 661 660 &namdrg_top ! TOP friction (ln_OFF =F & ln_isfcav=T) 662 661 !----------------------------------------------------------------------- 663 rn_Cd0 = 1.e-3! drag coefficient [-]664 rn_Uc0 = 0.4! ref. velocity [m/s] (linear drag=Cd0*Uc0)665 rn_Cdmax = 0.1! drag value maximum [-] (logarithmic drag)666 rn_ke0 = 2.5e-3! background kinetic energy [m2/s2] (non-linear cases)667 rn_z0 = 3.0e-3! roughness [m] (ln_loglayer=T)668 ln_boost = .false.! =T regional boost of Cd0 ; =F constant669 rn_boost = 50.! local boost factor [-]662 rn_Cd0 = 1.e-3 ! drag coefficient [-] 663 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 664 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag) 665 rn_ke0 = 2.5e-3 ! background kinetic energy [m2/s2] (non-linear cases) 666 rn_z0 = 3.0e-3 ! roughness [m] (ln_loglayer=T) 667 ln_boost = .false. ! =T regional boost of Cd0 ; =F constant 668 rn_boost = 50. ! local boost factor [-] 670 669 / 671 670 !----------------------------------------------------------------------- 672 671 &namdrg_bot ! BOTTOM friction (ln_OFF =F) 673 672 !----------------------------------------------------------------------- 674 rn_Cd0 = 1.e-3 ! drag coefficient [-]675 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0)676 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag)677 rn_ke0 = 2.5e-3 ! background kinetic energy [m2/s2] (non-linear cases)678 rn_z0 = 3.e-3 ! roughness [m] (ln_loglayer=T)679 ln_boost = .false. ! =T regional boost of Cd0 ; =F constant680 rn_boost = 50. ! local boost factor [-]673 rn_Cd0 = 1.e-3 ! drag coefficient [-] 674 rn_Uc0 = 0.4 ! ref. velocity [m/s] (linear drag=Cd0*Uc0) 675 rn_Cdmax = 0.1 ! drag value maximum [-] (logarithmic drag) 676 rn_ke0 = 2.5e-3 ! background kinetic energy [m2/s2] (non-linear cases) 677 rn_z0 = 3.e-3 ! roughness [m] (ln_loglayer=T) 678 ln_boost = .false. ! =T regional boost of Cd0 ; =F constant 679 rn_boost = 50. ! local boost factor [-] 681 680 / 682 681 !----------------------------------------------------------------------- … … 684 683 !----------------------------------------------------------------------- 685 684 ln_trabbc = .false. ! Apply a geothermal heating at the ocean bottom 686 nn_geoflx = 2 ! geothermal heat flux: = 1 constant flux687 ! ! = 2 read variable flux [mW/m2]688 rn_geoflx_cst = 86.4e-3 ! Constant value of geothermal heat flux [mW/m2]685 nn_geoflx = 2 ! geothermal heat flux: = 1 constant flux 686 ! ! = 2 read variable flux [mW/m2] 687 rn_geoflx_cst = 86.4e-3 ! Constant value of geothermal heat flux [mW/m2] 689 688 690 689 cn_dir = './' ! root directory for the geothermal data location … … 866 865 ln_dynvor_ens = .false. ! enstrophy conserving scheme 867 866 ln_dynvor_mix = .false. ! mixed scheme 867 ln_dynvor_enT = .false. ! energy conserving scheme (T-point) 868 ln_dynvor_eeT = .false. ! energy conserving scheme (een using e3t) 868 869 ln_dynvor_een = .false. ! energy & enstrophy scheme 869 nn_een_e3f = 1 ! =0 e3f = mean masked e3t divided by 4 870 ! ! =1 e3f = mean masked e3t divided by the sum of mask 871 ln_dynvor_msk = .false. ! vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) ! PLEASE DO NOT ACTIVATE 870 nn_een_e3f = 1 ! =0 e3f = mi(mj(e3t))/4 871 ! ! =1 e3f = mi(mj(e3t))/mi(mj( tmask)) 872 ln_dynvor_msk = .false. ! vorticity multiplied by fmask (=T) ==>>> PLEASE DO NOT ACTIVATE 873 ! ! (f-point vorticity schemes only) 872 874 / 873 875 !----------------------------------------------------------------------- -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r9506 r9528 36 36 USE sbcapr ! surface boundary condition: atmospheric pressure 37 37 USE dynadv , ONLY: ln_dynadv_vec 38 USE dynvor ! vortivity scheme indicators 38 39 USE phycst ! physical constants 39 40 USE dynvor ! vorticity term … … 103 104 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 104 105 ! 105 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 106 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 106 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 107 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 108 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 107 109 ! 108 110 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) … … 149 151 INTEGER :: ikbu, iktu, noffset ! local integers 150 152 INTEGER :: ikbv, iktv ! - - 151 REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars152 REAL(wp) :: zx1, zx2, zu_spg, zhura ! - -153 REAL(wp) :: zy1, zy2, zv_spg, zhvra ! - -154 REAL(wp) :: za0, za1, za2, za3 ! - -155 REAL(wp) :: zmdi, zztmp ! - -153 REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars 154 REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - - 155 REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - - 156 REAL(wp) :: za0, za1, za2, za3 ! - - 157 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 156 158 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 157 159 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 158 160 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 159 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e 161 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 160 162 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 161 163 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points … … 237 239 ! 238 240 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 239 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 241 ! 242 SELECT CASE( nvor_scheme ) 243 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme) 240 244 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 241 245 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) … … 250 254 DO jj = 1, jpjm1 251 255 DO ji = 1, jpim1 252 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 253 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) & 254 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 255 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 256 !!gm bug for ocean cavities 257 ! zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 258 ! & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) & 259 ! & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 260 ! & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 261 !! replace by : 262 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 263 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 264 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 265 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 266 !!gm end 256 267 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 257 268 END DO … … 270 281 END DO 271 282 ! 272 ELSE !== all other schemes (ENE, ENS, MIX) 283 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 284 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 285 DO jj = 2, jpj 286 DO ji = 2, jpi 287 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 288 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 289 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 290 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 291 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 292 END DO 293 END DO 294 ! 295 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! 296 ! 273 297 zwz(:,:) = 0._wp 274 298 zhf(:,:) = 0._wp … … 280 304 !!gm 281 305 !! 282 306 IF( .NOT.ln_sco ) THEN 283 307 284 308 !!gm agree the JC comment : this should be done in a much clear way … … 290 314 ! ENDIF 291 315 ! zhf(:,:) = gdepw_0(:,:,jk+1) 292 ELSE 293 !zhf(:,:) = hbatf(:,:) 294 DO jj = 1, jpjm1 295 DO ji = 1, jpim1 296 zhf(ji,jj) = MAX( 0._wp, & 297 & ( ht_0(ji ,jj )*tmask(ji ,jj ,1) + & 298 & ht_0(ji+1,jj )*tmask(ji+1,jj ,1) + & 299 & ht_0(ji ,jj+1)*tmask(ji ,jj+1,1) + & 300 & ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) / & 301 & ( tmask(ji ,jj ,1) + tmask(ji+1,jj ,1) +& 302 & tmask(ji ,jj+1,1) + tmask(ji+1,jj+1,1) +& 303 & rsmall ) ) 304 END DO 305 END DO 306 END IF 307 308 DO jj = 1, jpjm1 309 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 310 END DO 316 ! 317 ELSE 318 ! 319 !zhf(:,:) = hbatf(:,:) 320 DO jj = 1, jpjm1 321 DO ji = 1, jpim1 322 !!gm Bug here in case of ocean cavities, further more ht_0 is masked by definition ==>> remove mask multiplication 323 ! zhf(ji,jj) = MAX( 0._wp, & 324 ! & ( ht_0(ji ,jj )*tmask(ji ,jj ,1) + & 325 ! & ht_0(ji+1,jj )*tmask(ji+1,jj ,1) + & 326 ! & ht_0(ji ,jj+1)*tmask(ji ,jj+1,1) + & 327 ! & ht_0(ji+1,jj+1)*tmask(ji+1,jj+1,1) ) / & 328 ! & ( tmask(ji ,jj ,1) + tmask(ji+1,jj ,1) +& 329 ! & tmask(ji ,jj+1,1) + tmask(ji+1,jj+1,1) +& 330 ! & rsmall ) ) 331 !!gm replaced by : 332 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 333 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 334 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 335 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 311 336 !!gm end 312 337 END DO 338 END DO 339 ENDIF 340 ! 341 DO jj = 1, jpjm1 342 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 343 END DO 344 ! 313 345 DO jk = 1, jpkm1 314 346 DO jj = 1, jpjm1 … … 324 356 END DO 325 357 zwz(:,:) = ff_f(:,:) * zwz(:,:) 326 END IF358 END SELECT 327 359 ENDIF 328 360 ! … … 369 401 ! !* barotropic Coriolis trends (vorticity scheme dependent) 370 402 ! ! -------------------------------------------------------- 403 ! 371 404 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 372 405 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 373 406 ! 374 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme 407 SELECT CASE( nvor_scheme ) 408 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 409 DO jj = 2, jpjm1 410 DO ji = 2, jpim1 ! vector opt. 411 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj) & 412 & * ( e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) & 413 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 414 ! 415 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj) & 416 & * ( e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) & 417 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 418 END DO 419 END DO 420 ! 421 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 375 422 DO jj = 2, jpjm1 376 423 DO ji = fs_2, fs_jpim1 ! vector opt. … … 385 432 END DO 386 433 ! 387 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme434 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 388 435 DO jj = 2, jpjm1 389 436 DO ji = fs_2, fs_jpim1 ! vector opt. … … 397 444 END DO 398 445 ! 399 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme446 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 400 447 DO jj = 2, jpjm1 401 448 DO ji = fs_2, fs_jpim1 ! vector opt. … … 411 458 END DO 412 459 ! 413 END IF460 END SELECT 414 461 ! 415 462 ! !* Right-Hand-Side of the barotropic momentum equation 416 463 ! ! ---------------------------------------------------- 417 464 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 418 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 419 DO jj = 2, jpjm1 420 DO ji = 2, jpim1 421 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 422 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 423 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 465 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 466 DO jj = 2, jpjm1 467 DO ji = 2, jpim1 468 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 469 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 470 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 471 & > rn_wdmin1 + rn_wdmin2 472 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 473 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 474 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 475 IF(ll_tmp1) THEN 476 zcpx(ji,jj) = 1.0_wp 477 ELSEIF(ll_tmp2) THEN 478 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 479 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 480 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 481 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 482 ELSE 483 zcpx(ji,jj) = 0._wp 484 ENDIF 485 ! 486 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 487 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 488 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 424 489 & > rn_wdmin1 + rn_wdmin2 425 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 426 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 427 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 428 429 IF(ll_tmp1) THEN 430 zcpx(ji,jj) = 1.0_wp 431 ELSE IF(ll_tmp2) THEN 432 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 433 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 434 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 435 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 436 437 ELSE 438 zcpx(ji,jj) = 0._wp 439 END IF 440 441 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 442 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 443 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 444 & > rn_wdmin1 + rn_wdmin2 445 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 446 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 447 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 490 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 491 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 492 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 448 493 449 IF(ll_tmp1) THEN 450 zcpy(ji,jj) = 1.0_wp 451 ELSE IF(ll_tmp2) THEN 452 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 453 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 454 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 455 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 456 457 ELSE 458 zcpy(ji,jj) = 0._wp 459 END IF 460 END DO 461 END DO 462 463 DO jj = 2, jpjm1 464 DO ji = 2, jpim1 465 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 466 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 467 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 468 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 469 470 END DO 471 END DO 472 494 IF(ll_tmp1) THEN 495 zcpy(ji,jj) = 1.0_wp 496 ELSE IF(ll_tmp2) THEN 497 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 498 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 499 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 500 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 501 ELSE 502 zcpy(ji,jj) = 0._wp 503 ENDIF 504 END DO 505 END DO 506 ! 507 DO jj = 2, jpjm1 508 DO ji = 2, jpim1 509 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 510 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 511 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 512 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 513 END DO 514 END DO 515 ! 473 516 ELSE 474 517 ! … … 675 718 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 676 719 vn_adv(:,:) = 0._wp 677 ! ! ==================== ! 678 679 IF (ln_wd_dl) THEN 720 ! 721 IF( ln_wd_dl ) THEN 680 722 zuwdmask(:,:) = 0._wp ! set to zero for definiteness (not sure this is necessary) 681 723 zvwdmask(:,:) = 0._wp ! 682 zuwdav2 (:,:) =0._wp683 zvwdav2 (:,:) =0._wp724 zuwdav2 (:,:) = 0._wp 725 zvwdav2 (:,:) = 0._wp 684 726 END IF 685 727 686 728 ! ! ==================== ! 687 729 DO jn = 1, icycle ! sub-time-step loop ! 688 730 ! ! ==================== ! … … 715 757 ! set wetting & drying mask at tracer points for this barotropic sub-step 716 758 IF ( ln_wd_dl ) THEN 717 759 ! 718 760 IF ( ln_wd_dl_rmp ) THEN 719 761 DO jj = 1, jpj … … 736 778 ELSE 737 779 ztwdmask(ji,jj) = 0._wp 738 END 780 ENDIF 739 781 END DO 740 782 END DO 741 END IF 742 743 END IF 744 745 783 ENDIF 784 ! 785 ENDIF 786 ! 746 787 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 747 788 DO ji = 2, fs_jpim1 ! Vector opt. … … 756 797 CALL lbc_lnk_multi( zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 757 798 ! 758 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 759 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 799 zhup2_e(:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 800 zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 801 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 760 802 ELSE 761 zhup2_e (:,:) = hu_n(:,:) 762 zhvp2_e (:,:) = hv_n(:,:) 803 zhup2_e(:,:) = hu_n(:,:) 804 zhvp2_e(:,:) = hv_n(:,:) 805 zhtp2_e(:,:) = ht_n(:,:) 763 806 ENDIF 764 807 ! !* after ssh … … 795 838 ENDIF 796 839 #endif 797 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt)840 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 798 841 799 842 IF ( ln_wd_dl ) THEN 800 801 802 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 803 843 ! 844 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 845 ! 804 846 DO jj = 1, jpjm1 805 847 DO ji = 1, jpim1 … … 821 863 END DO 822 864 END DO 823 824 825 END IF 865 ! 866 ENDIF 826 867 827 868 ! Sum over sub-time-steps to compute advective velocities … … 896 937 ENDIF 897 938 ! 898 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 899 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 939 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 940 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 941 900 942 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 901 943 DO jj = 2, jpjm1 … … 917 959 ELSE 918 960 zcpx(ji,jj) = 0._wp 919 END 920 961 ENDIF 962 ! 921 963 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 922 964 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 929 971 IF(ll_tmp1) THEN 930 972 zcpy(ji,jj) = 1.0_wp 931 ELSE 973 ELSEIF(ll_tmp2) THEN 932 974 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 933 975 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & … … 935 977 ELSE 936 978 zcpy(ji,jj) = 0._wp 937 END 979 ENDIF 938 980 END DO 939 981 END DO 940 END 982 ENDIF 941 983 ! 942 984 ! Compute associated depths at U and V points: … … 955 997 END DO 956 998 END DO 957 999 ! 958 1000 ENDIF 959 1001 ! … … 965 1007 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 966 1008 ! 967 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 1009 SELECT CASE( nvor_scheme ) 1010 CASE( np_ENT ) ! energy conserving scheme (t-point) 1011 DO jj = 2, jpjm1 1012 DO ji = 2, jpim1 ! vector opt. 1013 1014 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 1015 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1016 1017 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1018 & * ( e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) ) & 1019 & + e1e2t(ji ,jj)*zhtp2_e(ji ,jj)*ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) ) 1020 ! 1021 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1022 & * ( e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) ) & 1023 & + e1e2t(ji,jj )*zhtp2_e(ji,jj )*ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) ) 1024 END DO 1025 END DO 1026 ! 1027 CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point) 968 1028 DO jj = 2, jpjm1 969 1029 DO ji = fs_2, fs_jpim1 ! vector opt. … … 977 1037 END DO 978 1038 ! 979 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!1039 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 980 1040 DO jj = 2, jpjm1 981 1041 DO ji = fs_2, fs_jpim1 ! vector opt. … … 989 1049 END DO 990 1050 ! 991 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!1051 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 992 1052 DO jj = 2, jpjm1 993 1053 DO ji = fs_2, fs_jpim1 ! vector opt. 994 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &995 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &996 & + ftse(ji,jj ) * zwy(ji ,jj-1) &997 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )998 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &999 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &1000 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &1001 & + ftne(ji,jj ) * zwx(ji ,jj ) )1054 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 1055 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 1056 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 1057 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 1058 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 1059 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 1060 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 1061 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 1002 1062 END DO 1003 1063 END DO 1004 1064 ! 1005 END IF1065 END SELECT 1006 1066 ! 1007 1067 ! Add tidal astronomical forcing if defined … … 1341 1401 !! ** Purpose : Read or write time-splitting arrays in restart file 1342 1402 !!---------------------------------------------------------------------- 1343 INTEGER , INTENT(in) :: kt ! ocean time-step 1344 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1345 ! 1346 INTEGER :: id1, id2 ! local integers 1403 INTEGER , INTENT(in) :: kt ! ocean time-step 1404 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1347 1405 !!---------------------------------------------------------------------- 1348 1406 ! -
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r9190 r9528 6 6 !!====================================================================== 7 7 !! History : OPA ! 1989-12 (P. Andrich) vor_ens: Original code 8 !! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code8 !! 5.0 ! 1991-11 (G. Madec) vor_ene, vor_mix: Original code 9 9 !! 6.0 ! 1996-01 (G. Madec) s-coord, suppress work arrays 10 10 !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module … … 19 19 !! - ! 2016-12 (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) 20 20 !! 4.0 ! 2017-07 (G. Madec) linear dynamics + trends diag. with Stokes-Coriolis 21 !! - ! 2018-03 (G. Madec) add two new schemes (ln_dynvor_enT and ln_dynvor_eet) 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 21 23 !!---------------------------------------------------------------------- 22 24 … … 50 52 51 53 ! !!* Namelist namdyn_vor: vorticity term 52 LOGICAL, PUBLIC :: ln_dynvor_ene !: energy conserving scheme (ENE) 53 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 54 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 55 LOGICAL, PUBLIC :: ln_dynvor_een !: energy and enstrophy conserving scheme (EEN) 54 LOGICAL, PUBLIC :: ln_dynvor_ens !: enstrophy conserving scheme (ENS) 55 LOGICAL, PUBLIC :: ln_dynvor_ene !: f-point energy conserving scheme (ENE) 56 LOGICAL, PUBLIC :: ln_dynvor_enT !: t-point energy conserving scheme (ENT) 57 LOGICAL, PUBLIC :: ln_dynvor_eeT !: t-point energy conserving scheme (EET) 58 LOGICAL, PUBLIC :: ln_dynvor_een !: energy & enstrophy conserving scheme (EEN) 56 59 INTEGER, PUBLIC :: nn_een_e3f !: e3f=masked averaging of e3t divided by 4 (=0) or by the sum of mask (=1) 60 LOGICAL, PUBLIC :: ln_dynvor_mix !: mixed scheme (MIX) 57 61 LOGICAL, PUBLIC :: ln_dynvor_msk !: vorticity multiplied by fmask (=T) or not (=F) (all vorticity schemes) 58 62 59 INTEGER :: nvor_scheme ! choice of the type of advection scheme 60 ! ! associated indices: 63 INTEGER, PUBLIC :: nvor_scheme !: choice of the type of advection scheme 64 ! ! associated indices: 65 INTEGER, PUBLIC, PARAMETER :: np_ENS = 0 ! ENS scheme 61 66 INTEGER, PUBLIC, PARAMETER :: np_ENE = 1 ! ENE scheme 62 INTEGER, PUBLIC, PARAMETER :: np_EN S = 2 ! ENS scheme63 INTEGER, PUBLIC, PARAMETER :: np_ MIX = 3 ! MIX scheme67 INTEGER, PUBLIC, PARAMETER :: np_ENT = 2 ! ENT scheme (t-point vorticity) 68 INTEGER, PUBLIC, PARAMETER :: np_EET = 3 ! EET scheme (EEN using e3t) 64 69 INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme 70 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 65 71 66 72 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 67 73 ! ! associated indices: 68 INTEGER, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 69 INTEGER, PARAMETER :: np_RVO = 2 ! relative vorticity 70 INTEGER, PARAMETER :: np_MET = 3 ! metric term 71 INTEGER, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 72 INTEGER, PARAMETER :: np_CME = 5 ! Coriolis + metric term 74 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 75 INTEGER, PUBLIC, PARAMETER :: np_RVO = 2 ! relative vorticity 76 INTEGER, PUBLIC, PARAMETER :: np_MET = 3 ! metric term 77 INTEGER, PUBLIC, PARAMETER :: np_CRV = 4 ! relative + planetary (total vorticity) 78 INTEGER, PUBLIC, PARAMETER :: np_CME = 5 ! Coriolis + metric term 79 80 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2u_2 ! = di(e2u)/2 used in T-point metric term calculation 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - - 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: di_e2v_2e1e2f ! = di(e2u)/(2*e1e2f) used in F-point metric term calculation 83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1v)/(2*e1e2f) - - - - 73 84 74 85 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 … … 109 120 ztrdv(:,:,:) = va(:,:,:) 110 121 SELECT CASE( nvor_scheme ) 122 CASE( np_ENS ) ; CALL vor_ens( kt, ncor, un , vn , ua, va ) ! enstrophy conserving scheme 123 IF( ln_stcor ) CALL vor_ens( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 111 124 CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, ncor, un , vn , ua, va ) ! energy conserving scheme 112 125 IF( ln_stcor ) CALL vor_ene( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 113 CASE( np_ENS ) ; CALL vor_ens( kt, ncor, un , vn , ua, va ) ! enstrophy conserving scheme 114 IF( ln_stcor ) CALL vor_ens( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 126 CASE( np_ENT ) ; CALL vor_enT( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (T-pts) 127 IF( ln_stcor ) CALL vor_enT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 128 CASE( np_EET ) ; CALL vor_eeT( kt, ncor, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 129 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 115 130 CASE( np_EEN ) ; CALL vor_een( kt, ncor, un , vn , ua, va ) ! energy & enstrophy scheme 116 131 IF( ln_stcor ) CALL vor_een( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend … … 124 139 ztrdv(:,:,:) = va(:,:,:) 125 140 SELECT CASE( nvor_scheme ) 141 CASE( np_ENT ) ; CALL vor_enT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (T-pts) 142 CASE( np_EET ) ; CALL vor_eeT( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme (een with e3t) 126 143 CASE( np_ENE ) ; CALL vor_ene( kt, nrvm, un , vn , ua, va ) ! energy conserving scheme 127 144 CASE( np_ENS, np_MIX ) ; CALL vor_ens( kt, nrvm, un , vn , ua, va ) ! enstrophy conserving scheme … … 138 155 ! 139 156 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 157 CASE( np_ENT ) !* energy conserving scheme (T-pts) 158 CALL vor_enT( kt, ntot, un , vn , ua, va ) ! total vorticity trend 159 IF( ln_stcor ) CALL vor_enT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 160 CASE( np_EET ) !* energy conserving scheme (een scheme using e3t) 161 CALL vor_eeT( kt, ntot, un , vn , ua, va ) ! total vorticity trend 162 IF( ln_stcor ) CALL vor_eeT( kt, ncor, usd, vsd, ua, va ) ! add the Stokes-Coriolis trend 140 163 CASE( np_ENE ) !* energy conserving scheme 141 164 CALL vor_ene( kt, ntot, un , vn , ua, va ) ! total vorticity trend … … 164 187 165 188 189 SUBROUTINE vor_enT( kt, kvor, pu, pv, pu_rhs, pv_rhs ) 190 !!---------------------------------------------------------------------- 191 !! *** ROUTINE vor_enT *** 192 !! 193 !! ** Purpose : Compute the now total vorticity trend and add it to 194 !! the general trend of the momentum equation. 195 !! 196 !! ** Method : Trend evaluated using now fields (centered in time) 197 !! and t-point evaluation of vorticity (planetary and relative). 198 !! conserves the horizontal kinetic energy. 199 !! The general trend of momentum is increased due to the vorticity 200 !! term which is given by: 201 !! voru = 1/bu mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t mj[vn] ] 202 !! vorv = 1/bv mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f mj[un] ] 203 !! where rvor is the relative vorticity at f-point 204 !! 205 !! ** Action : - Update (ua,va) with the now vorticity term trend 206 !!---------------------------------------------------------------------- 207 INTEGER , INTENT(in ) :: kt ! ocean time-step index 208 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 209 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities 210 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend 211 ! 212 INTEGER :: ji, jj, jk ! dummy loop indices 213 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 214 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zwt ! 2D workspace 215 !!---------------------------------------------------------------------- 216 ! 217 IF( kt == nit000 ) THEN 218 IF(lwp) WRITE(numout,*) 219 IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme' 220 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 221 ENDIF 222 ! 223 ! ! =============== 224 DO jk = 1, jpkm1 ! Horizontal slab 225 ! ! =============== 226 ! 227 SELECT CASE( kvor ) !== volume weighted vorticity considered ==! 228 CASE ( np_COR ) !* Coriolis (planetary vorticity) 229 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t_n(:,:,jk) 230 CASE ( np_RVO ) !* relative vorticity 231 DO jj = 1, jpjm1 232 DO ji = 1, jpim1 233 zwz(ji,jj) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 234 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 235 END DO 236 END DO 237 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 238 DO jj = 1, jpjm1 239 DO ji = 1, jpim1 240 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 241 END DO 242 END DO 243 ENDIF 244 CALL lbc_lnk( zwz, 'F', 1. ) 245 DO jj = 2, jpj 246 DO ji = 2, jpi ! vector opt. 247 zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) & 248 & + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 249 END DO 250 END DO 251 CASE ( np_MET ) !* metric term 252 DO jj = 2, jpj 253 DO ji = 2, jpi 254 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 255 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t_n(ji,jj,jk) 256 END DO 257 END DO 258 CASE ( np_CRV ) !* Coriolis + relative vorticity 259 DO jj = 1, jpjm1 260 DO ji = 1, jpim1 ! relative vorticity 261 zwz(ji,jj) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 262 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 263 END DO 264 END DO 265 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 266 DO jj = 1, jpjm1 267 DO ji = 1, jpim1 268 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 269 END DO 270 END DO 271 ENDIF 272 CALL lbc_lnk( zwz, 'F', 1. ) 273 DO jj = 2, jpj 274 DO ji = 2, jpi ! vector opt. 275 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) & 276 & + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) ) * e1e2t(ji,jj)*e3t_n(ji,jj,jk) 277 END DO 278 END DO 279 CASE ( np_CME ) !* Coriolis + metric 280 DO jj = 2, jpj 281 DO ji = 2, jpi ! vector opt. 282 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 283 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 284 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t_n(ji,jj,jk) 285 END DO 286 END DO 287 CASE DEFAULT ! error 288 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 289 END SELECT 290 ! 291 ! !== compute and add the vorticity term trend =! 292 DO jj = 2, jpjm1 293 DO ji = 2, jpim1 ! vector opt. 294 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) & 295 & * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) & 296 & + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) ) 297 ! 298 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) & 299 & * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) & 300 & + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) ) 301 END DO 302 END DO 303 ! ! =============== 304 END DO ! End of slab 305 ! ! =============== 306 END SUBROUTINE vor_enT 307 308 166 309 SUBROUTINE vor_ene( kt, kvor, pun, pvn, pua, pva ) 167 310 !!---------------------------------------------------------------------- … … 217 360 DO jj = 1, jpjm1 218 361 DO ji = 1, fs_jpim1 ! vector opt. 219 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 220 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 221 & * 0.5 * r1_e1e2f(ji,jj) 362 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 363 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 222 364 END DO 223 365 END DO … … 225 367 DO jj = 1, jpjm1 226 368 DO ji = 1, fs_jpim1 ! vector opt. 227 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 228 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 229 & * r1_e1e2f(ji,jj) 369 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 370 & - e1u(ji,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 230 371 END DO 231 372 END DO … … 233 374 DO jj = 1, jpjm1 234 375 DO ji = 1, fs_jpim1 ! vector opt. 235 zwz(ji,jj) = ff_f(ji,jj) & 236 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 237 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 238 & * 0.5 * r1_e1e2f(ji,jj) 376 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 377 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 239 378 END DO 240 379 END DO … … 328 467 DO jj = 1, jpjm1 329 468 DO ji = 1, fs_jpim1 ! vector opt. 330 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 331 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 332 & * 0.5 * r1_e1e2f(ji,jj) 469 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 470 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 333 471 END DO 334 472 END DO … … 336 474 DO jj = 1, jpjm1 337 475 DO ji = 1, fs_jpim1 ! vector opt. 338 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 339 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 340 & * r1_e1e2f(ji,jj) 476 zwz(ji,jj) = ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 477 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj) 341 478 END DO 342 479 END DO … … 344 481 DO jj = 1, jpjm1 345 482 DO ji = 1, fs_jpim1 ! vector opt. 346 zwz(ji,jj) = ff_f(ji,jj) & 347 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 348 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 349 & * 0.5 * r1_e1e2f(ji,jj) 483 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 484 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 350 485 END DO 351 486 END DO … … 412 547 INTEGER :: ierr ! local integer 413 548 REAL(wp) :: zua, zva ! local scalars 414 REAL(wp) :: zmsk, ze3 415 REAL(wp), DIMENSION(jpi,jpj) ::zwx , zwy , zwz , z1_e3f416 REAL(wp), DIMENSION(jpi,jpj) ::ztnw, ztne, ztsw, ztse549 REAL(wp) :: zmsk, ze3f ! local scalars 550 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , zwz , z1_e3f 551 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 417 552 !!---------------------------------------------------------------------- 418 553 ! … … 431 566 DO jj = 1, jpjm1 432 567 DO ji = 1, fs_jpim1 ! vector opt. 433 ze3 568 ze3f = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 434 569 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 435 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3436 ELSE ; z1_e3f(ji,jj) = 0._wp570 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f 571 ELSE ; z1_e3f(ji,jj) = 0._wp 437 572 ENDIF 438 573 END DO … … 441 576 DO jj = 1, jpjm1 442 577 DO ji = 1, fs_jpim1 ! vector opt. 443 ze3 578 ze3f = ( e3t_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk) & 444 579 & + e3t_n(ji,jj ,jk)*tmask(ji,jj ,jk) + e3t_n(ji+1,jj ,jk)*tmask(ji+1,jj ,jk) ) 445 580 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 446 581 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) 447 IF( ze3 /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3448 ELSE ; z1_e3f(ji,jj) = 0._wp582 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f 583 ELSE ; z1_e3f(ji,jj) = 0._wp 449 584 ENDIF 450 585 END DO … … 462 597 DO jj = 1, jpjm1 463 598 DO ji = 1, fs_jpim1 ! vector opt. 464 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 465 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 466 & * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 599 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 600 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 467 601 END DO 468 602 END DO … … 470 604 DO jj = 1, jpjm1 471 605 DO ji = 1, fs_jpim1 ! vector opt. 472 zwz(ji,jj) = ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 473 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 474 & * 0.5 * r1_e1e2f(ji,jj) * z1_e3f(ji,jj) 606 zwz(ji,jj) = ( ( pvn(ji+1,jj,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 607 & - ( pun(ji,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 475 608 END DO 476 609 END DO … … 478 611 DO jj = 1, jpjm1 479 612 DO ji = 1, fs_jpim1 ! vector opt. 480 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk)&481 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) )&482 & * r1_e1e2f(ji,jj)) * z1_e3f(ji,jj)613 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 614 & - e1u(ji ,jj+1) * pun(ji,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 615 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 483 616 END DO 484 617 END DO … … 486 619 DO jj = 1, jpjm1 487 620 DO ji = 1, fs_jpim1 ! vector opt. 488 zwz(ji,jj) = ( ff_f(ji,jj) & 489 & + ( ( pvn(ji+1,jj ,jk) + pvn (ji,jj,jk) ) * ( e2v(ji+1,jj ) - e2v(ji,jj) ) & 490 & - ( pun(ji ,jj+1,jk) + pun (ji,jj,jk) ) * ( e1u(ji ,jj+1) - e1u(ji,jj) ) ) & 491 & * 0.5 * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 621 zwz(ji,jj) = ( ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 622 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 492 623 END DO 493 624 END DO … … 543 674 544 675 676 677 SUBROUTINE vor_eeT( kt, kvor, pun, pvn, pua, pva ) 678 !!---------------------------------------------------------------------- 679 !! *** ROUTINE vor_eeT *** 680 !! 681 !! ** Purpose : Compute the now total vorticity trend and add it to 682 !! the general trend of the momentum equation. 683 !! 684 !! ** Method : Trend evaluated using now fields (centered in time) 685 !! and the Arakawa and Lamb (1980) vector form formulation using 686 !! a modified version of Arakawa and Lamb (1980) scheme (see vor_een). 687 !! The change consists in 688 !! Add this trend to the general momentum trend (ua,va). 689 !! 690 !! ** Action : - Update (ua,va) with the now vorticity term trend 691 !! 692 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 693 !!---------------------------------------------------------------------- 694 INTEGER , INTENT(in ) :: kt ! ocean time-step index 695 INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric 696 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pun, pvn ! now velocities 697 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pua, pva ! total v-trend 698 ! 699 INTEGER :: ji, jj, jk ! dummy loop indices 700 INTEGER :: ierr ! local integer 701 REAL(wp) :: zua, zva ! local scalars 702 REAL(wp) :: zmsk, z1_e3t ! local scalars 703 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , zwz 704 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 705 !!---------------------------------------------------------------------- 706 ! 707 IF( kt == nit000 ) THEN 708 IF(lwp) WRITE(numout,*) 709 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 710 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 711 ENDIF 712 ! 713 ! ! =============== 714 DO jk = 1, jpkm1 ! Horizontal slab 715 ! ! =============== 716 ! 717 ! 718 SELECT CASE( kvor ) !== vorticity considered ==! 719 CASE ( np_COR ) !* Coriolis (planetary vorticity) 720 DO jj = 1, jpjm1 721 DO ji = 1, fs_jpim1 ! vector opt. 722 zwz(ji,jj) = ff_f(ji,jj) 723 END DO 724 END DO 725 CASE ( np_RVO ) !* relative vorticity 726 DO jj = 1, jpjm1 727 DO ji = 1, fs_jpim1 ! vector opt. 728 zwz(ji,jj) = ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 729 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 730 & * r1_e1e2f(ji,jj) 731 END DO 732 END DO 733 CASE ( np_MET ) !* metric term 734 DO jj = 1, jpjm1 735 DO ji = 1, fs_jpim1 ! vector opt. 736 zwz(ji,jj) = ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 737 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 738 END DO 739 END DO 740 CASE ( np_CRV ) !* Coriolis + relative vorticity 741 DO jj = 1, jpjm1 742 DO ji = 1, fs_jpim1 ! vector opt. 743 zwz(ji,jj) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pvn(ji+1,jj ,jk) - e2v(ji,jj) * pvn(ji,jj,jk) & 744 & - e1u(ji ,jj+1) * pun(ji ,jj+1,jk) + e1u(ji,jj) * pun(ji,jj,jk) ) & 745 & * r1_e1e2f(ji,jj) ) 746 END DO 747 END DO 748 CASE ( np_CME ) !* Coriolis + metric 749 DO jj = 1, jpjm1 750 DO ji = 1, fs_jpim1 ! vector opt. 751 zwz(ji,jj) = ff_f(ji,jj) + ( pvn(ji+1,jj ,jk) + pvn(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 752 & - ( pun(ji ,jj+1,jk) + pun(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 753 END DO 754 END DO 755 CASE DEFAULT ! error 756 CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' ) 757 END SELECT 758 ! 759 IF( ln_dynvor_msk ) THEN !== mask/unmask vorticity ==! 760 DO jj = 1, jpjm1 761 DO ji = 1, fs_jpim1 ! vector opt. 762 zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk) 763 END DO 764 END DO 765 ENDIF 766 ! 767 CALL lbc_lnk( zwz, 'F', 1. ) 768 ! 769 ! !== horizontal fluxes ==! 770 zwx(:,:) = e2u(:,:) * e3u_n(:,:,jk) * pun(:,:,jk) 771 zwy(:,:) = e1v(:,:) * e3v_n(:,:,jk) * pvn(:,:,jk) 772 773 ! !== compute and add the vorticity term trend =! 774 jj = 2 775 ztne(1,:) = 0 ; ztnw(1,:) = 0 ; ztse(1,:) = 0 ; ztsw(1,:) = 0 776 DO ji = 2, jpi ! split in 2 parts due to vector opt. 777 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 778 ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) * z1_e3t 779 ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) * z1_e3t 780 ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t 781 ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) * z1_e3t 782 END DO 783 DO jj = 3, jpj 784 DO ji = fs_2, jpi ! vector opt. ok because we start at jj = 3 785 z1_e3t = 1._wp / e3t_n(ji,jj,jk) 786 ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) * z1_e3t 787 ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) * z1_e3t 788 ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t 789 ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) * z1_e3t 790 END DO 791 END DO 792 DO jj = 2, jpjm1 793 DO ji = fs_2, fs_jpim1 ! vector opt. 794 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & 795 & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) ) 796 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) & 797 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 798 pua(ji,jj,jk) = pua(ji,jj,jk) + zua 799 pva(ji,jj,jk) = pva(ji,jj,jk) + zva 800 END DO 801 END DO 802 ! ! =============== 803 END DO ! End of slab 804 ! ! =============== 805 END SUBROUTINE vor_eeT 806 807 545 808 SUBROUTINE dyn_vor_init 546 809 !!--------------------------------------------------------------------- … … 550 813 !! tracer advection schemes 551 814 !!---------------------------------------------------------------------- 552 INTEGER :: ioptio ! local integer 553 INTEGER :: ji, jj, jk ! dummy loop indices 554 INTEGER :: ios ! Local integer output status for namelist read 555 !! 556 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_mix, & 557 & ln_dynvor_een, nn_een_e3f , ln_dynvor_msk 558 !!---------------------------------------------------------------------- 559 815 INTEGER :: ji, jj, jk ! dummy loop indices 816 INTEGER :: ioptio, ios ! local integer 817 !! 818 NAMELIST/namdyn_vor/ ln_dynvor_ens, ln_dynvor_ene, ln_dynvor_enT, ln_dynvor_eeT, & 819 & ln_dynvor_een, nn_een_e3f , ln_dynvor_mix, ln_dynvor_msk 820 !!---------------------------------------------------------------------- 821 ! 822 IF(lwp) THEN 823 WRITE(numout,*) 824 WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency' 825 WRITE(numout,*) '~~~~~~~~~~~~' 826 ENDIF 827 ! 560 828 REWIND( numnam_ref ) ! Namelist namdyn_vor in reference namelist : Vorticity scheme options 561 829 READ ( numnam_ref, namdyn_vor, IOSTAT = ios, ERR = 901) … … 565 833 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_vor in configuration namelist', lwp ) 566 834 IF(lwm) WRITE ( numond, namdyn_vor ) 567 835 ! 568 836 IF(lwp) THEN ! Namelist print 569 WRITE(numout,*)570 WRITE(numout,*) 'dyn_vor_init : vorticity term : read namelist and control the consistency'571 WRITE(numout,*) '~~~~~~~~~~~~'572 837 WRITE(numout,*) ' Namelist namdyn_vor : choice of the vorticity term scheme' 573 WRITE(numout,*) ' energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene574 838 WRITE(numout,*) ' enstrophy conserving scheme ln_dynvor_ens = ', ln_dynvor_ens 575 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 839 WRITE(numout,*) ' f-point energy conserving scheme ln_dynvor_ene = ', ln_dynvor_ene 840 WRITE(numout,*) ' t-point energy conserving scheme ln_dynvor_enT = ', ln_dynvor_enT 841 WRITE(numout,*) ' energy conserving scheme (een using e3t) ln_dynvor_eeT = ', ln_dynvor_eeT 576 842 WRITE(numout,*) ' enstrophy and energy conserving scheme ln_dynvor_een = ', ln_dynvor_een 577 843 WRITE(numout,*) ' e3f = averaging /4 (=0) or /sum(tmask) (=1) nn_een_e3f = ', nn_een_e3f 844 WRITE(numout,*) ' mixed enstrophy/energy conserving scheme ln_dynvor_mix = ', ln_dynvor_mix 578 845 WRITE(numout,*) ' masked (=T) or unmasked(=F) vorticity ln_dynvor_msk = ', ln_dynvor_msk 579 846 ENDIF 847 848 IF( ln_dynvor_msk ) CALL ctl_stop( 'dyn_vor_init: masked vorticity is not currently not available') 580 849 581 850 !!gm this should be removed when choosing a unique strategy for fmask at the coast … … 586 855 IF( ln_vorlat .AND. ( ln_dynvor_ene .OR. ln_dynvor_ens .OR. ln_dynvor_mix ) ) THEN 587 856 DO jk = 1, jpk 588 DO jj = 2, jpjm1589 DO ji = 2, jpim1590 IF( tmask(ji,jj,jk)+tmask(ji+1,jj,jk)+tmask(ji,jj+1,jk)+tmask(ji+1,jj+1,jk) == 3._wp )&591 fmask(ji,jj,jk) = 1._wp857 DO jj = 1, jpjm1 858 DO ji = 1, jpim1 859 IF( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 860 & + tmask(ji,jj ,jk) + tmask(ji+1,jj+1,jk) == 3._wp ) fmask(ji,jj,jk) = 1._wp 592 861 END DO 593 862 END DO 594 863 END DO 595 596 597 864 ! 865 CALL lbc_lnk( fmask, 'F', 1._wp ) ! Lateral boundary conditions on fmask 866 ! 598 867 ENDIF 599 868 !!gm end 600 869 601 870 ioptio = 0 ! type of scheme for vorticity (set nvor_scheme) 602 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 603 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 604 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 605 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 871 IF( ln_dynvor_ens ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENS ; ENDIF 872 IF( ln_dynvor_ene ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENE ; ENDIF 873 IF( ln_dynvor_enT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_ENT ; ENDIF 874 IF( ln_dynvor_eeT ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EET ; ENDIF 875 IF( ln_dynvor_een ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_EEN ; ENDIF 876 IF( ln_dynvor_mix ) THEN ; ioptio = ioptio + 1 ; nvor_scheme = np_MIX ; ENDIF 606 877 ! 607 878 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) … … 622 893 nrvm = np_MET ! metric term 623 894 ntot = np_CME ! Coriolis + metric term 895 ! 896 SELECT CASE( nvor_scheme ) ! pre-computed gradients for the metric term: 897 CASE( np_ENT ) !* T-point metric term : pre-compute di(e2u)/2 and dj(e1v)/2 898 ALLOCATE( di_e2u_2(jpi,jpj), dj_e1v_2(jpi,jpj) ) 899 DO jj = 2, jpjm1 900 DO ji = 2, jpim1 901 di_e2u_2(ji,jj) = ( e2u(ji,jj) - e2u(ji-1,jj ) ) * 0.5_wp 902 dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp 903 END DO 904 END DO 905 CALL lbc_lnk_multi( di_e2u_2, 'T', -1. , dj_e1v_2, 'T', -1. ) ! Lateral boundary conditions 906 ! 907 CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f) 908 ALLOCATE( di_e2v_2e1e2f(jpi,jpj), dj_e1u_2e1e2f(jpi,jpj) ) 909 DO jj = 1, jpjm1 910 DO ji = 1, jpim1 911 di_e2v_2e1e2f(ji,jj) = ( e2v(ji+1,jj ) - e2v(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 912 dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 913 END DO 914 END DO 915 CALL lbc_lnk_multi( di_e2v_2e1e2f, 'F', -1. , dj_e1u_2e1e2f, 'F', -1. ) ! Lateral boundary conditions 916 END SELECT 917 ! 624 918 END SELECT 625 919 … … 627 921 WRITE(numout,*) 628 922 SELECT CASE( nvor_scheme ) 629 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme' 630 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme' 631 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme' 632 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme' 923 CASE( np_ENS ) ; WRITE(numout,*) ' ==>>> enstrophy conserving scheme (ENS)' 924 CASE( np_ENE ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at F-points) (ENE)' 925 CASE( np_ENT ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (Coriolis at T-points) (ENT)' 926 CASE( np_EET ) ; WRITE(numout,*) ' ==>>> energy conserving scheme (EEN scheme using e3t) (EET)' 927 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 928 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 633 929 END SELECT 634 930 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
