Changeset 14547
- Timestamp:
- 2021-02-25T18:07:15+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE
- Files:
-
- 3 added
- 11 edited
-
DOM/domain.F90 (modified) (2 diffs)
-
DYN/divhor.F90 (modified) (5 diffs)
-
DYN/dynspg.F90 (modified) (4 diffs)
-
DYN/dynspg_ts.F90 (modified) (31 diffs)
-
DYN/dynvor.F90 (modified) (5 diffs)
-
DYN/dynzdf.F90 (modified) (23 diffs)
-
DYN/sshwzv.F90 (modified) (9 diffs)
-
TRA/eosbn2.F90 (modified) (2 diffs)
-
TRA/tranpc.F90 (modified) (2 diffs)
-
nemogcm.F90 (modified) (3 diffs)
-
stp2d.F90 (added)
-
stpmlf.F90 (modified) (1 diff)
-
stprk3.F90 (added)
-
stprk3_stg.F90 (added)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DOM/domain.F90
r14255 r14547 91 91 INTEGER :: ji, jj, jk, jt ! dummy loop indices 92 92 INTEGER :: iconf = 0 ! local integers 93 REAL(wp):: zrdt94 93 CHARACTER (len=64) :: cform = "(A12, 3(A13, I7))" 95 94 INTEGER , DIMENSION(jpi,jpj) :: ik_top , ik_bot ! top and bottom ocean level … … 372 371 ! 373 372 ! set current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3 373 #if defined key_RK3 374 rDt = rn_Dt 375 r1_Dt = 1._wp / rDt 376 ! 377 IF(lwp) THEN 378 WRITE(numout,*) 379 WRITE(numout,*) ' ===>>> Runge Kutta 3rd order (RK3) : rDt = ', rDt 380 WRITE(numout,*) 381 ENDIF 382 ! 383 #else 374 384 rDt = 2._wp * rn_Dt 375 385 r1_Dt = 1._wp / rDt 386 ! 387 IF(lwp) THEN 388 WRITE(numout,*) 389 WRITE(numout,*) ' ===>>> Modified Leap-Frog (MLF) : rDt = ', rDt 390 WRITE(numout,*) 391 ENDIF 392 ! 393 #endif 376 394 ! 377 395 IF( l_SAS .AND. .NOT.ln_linssh ) THEN -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/divhor.F90
r13558 r14547 13 13 !! - ! 2014-12 (G. Madec) suppression of cross land advection option 14 14 !! - ! 2015-10 (G. Madec) add velocity and rnf flag in argument of div_hor 15 !! 4.5 ! 2015-10 (S. Techene, G. Madec) hdiv replaced by e3divUh 15 16 !!---------------------------------------------------------------------- 16 17 … … 36 37 PRIVATE 37 38 38 PUBLIC div_hor ! routine called by step.F90 and istate.F90 39 ! !! * Interface 40 INTERFACE div_hor 41 MODULE PROCEDURE div_hor_RK3, div_hor_old 42 END INTERFACE 43 44 PUBLIC div_hor ! routine called by ssh_nxt.F90 and istate.F90 39 45 40 46 !! * Substitutions … … 48 54 CONTAINS 49 55 50 SUBROUTINE div_hor( kt, Kbb, Kmm ) 56 SUBROUTINE div_hor_RK3( kt, Kbb, Kmm, puu, pvv, pe3divUh ) 57 !!---------------------------------------------------------------------- 58 !! *** ROUTINE div_hor_RK3 *** 59 !! 60 !! ** Purpose : compute the horizontal divergence at now time-step 61 !! 62 !! ** Method : the now divergence is computed as : 63 !! hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 64 !! and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 65 !! 66 !! ** Action : - thickness weighted horizontal divergence of in input velocity (puu,pvv) 67 !!---------------------------------------------------------------------- 68 INTEGER , INTENT(in ) :: kt, Kbb, Kmm ! ocean time-step & time-level indices 69 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: puu, pvv ! horizontal velocity 70 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: pe3divUh ! e3t*div[Uh] 71 ! 72 INTEGER :: ji, jj, jk ! dummy loop indices 73 !!---------------------------------------------------------------------- 74 ! 75 IF( ln_timing ) CALL timing_start('div_hor_RK3') 76 ! 77 IF( kt == nit000 ) THEN 78 IF(lwp) WRITE(numout,*) 79 IF(lwp) WRITE(numout,*) 'div_hor_RK3 : thickness weighted horizontal divergence ' 80 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 81 hdiv (:,:,:) = 0._wp ! initialize hdiv & pe3divUh for the halos and jpk level at the first time step 82 ENDIF 83 ! 84 pe3divUh(:,:,:) = 0._wp !!gm to be applied to the halos only 85 ! 86 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 87 hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * puu(ji ,jj,jk) & 88 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk) & 89 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * pvv(ji,jj ,jk) & 90 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * pvv(ji,jj-1,jk) ) & 91 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 92 END_3D 93 ! 94 IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== + runoffs divergence ==! 95 ! 96 #if defined key_asminc 97 IF( ln_sshinc .AND. ln_asmiau ) & !== + SSH assimilation increment ==! 98 & CALL ssh_asm_div( kt, Kbb, Kmm, hdiv ) 99 #endif 100 ! 101 IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== + ice-shelf mass exchange ==! 102 ! 103 !!gm Patch before suppression of hdiv from all modules that use it 104 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== e3t * Horizontal divergence ==! 105 pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 106 END_3D 107 !!gm end 108 ! 109 CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp ) ! (no sign change) 110 ! 111 IF( ln_timing ) CALL timing_stop('div_hor_RK3') 112 ! 113 END SUBROUTINE div_hor_RK3 114 115 116 SUBROUTINE div_hor_old( kt, Kbb, Kmm ) 51 117 !!---------------------------------------------------------------------- 52 118 !! *** ROUTINE div_hor *** … … 64 130 ! 65 131 INTEGER :: ji, jj, jk ! dummy loop indices 66 REAL(wp) :: zraur, zdep ! local scalars67 REAL(wp), DIMENSION(jpi,jpj) :: ztmp68 132 !!---------------------------------------------------------------------- 69 133 ! … … 98 162 IF( ln_timing ) CALL timing_stop('div_hor') 99 163 ! 100 END SUBROUTINE div_hor 164 END SUBROUTINE div_hor_old 101 165 102 166 !!====================================================================== -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg.F90
r14225 r14547 57 57 CONTAINS 58 58 59 SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa , k_only_ADV)59 SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 60 60 !!---------------------------------------------------------------------- 61 61 !! *** ROUTINE dyn_spg *** … … 79 79 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 80 80 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels 81 INTEGER , OPTIONAL , INTENT( in ) :: k_only_ADV ! only Advection in the RHS82 81 ! 83 82 INTEGER :: ji, jj, jk ! dummy loop indices 84 REAL(wp) :: z 2dt, zg_2, zintp, zgrho0r, zld! local scalars83 REAL(wp) :: zg_2, zintp, zgrho0r, zld ! local scalars 85 84 REAL(wp) , DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace 86 85 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice … … 151 150 IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head 152 151 DO_2D( 0, 0, 0, 0 ) 153 zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) /e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]154 zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) /e2v(ji,jj)152 zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3] 153 zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e2v(ji,jj) 155 154 END_2D 156 155 ENDIF … … 167 166 SELECT CASE ( nspg ) !== surface pressure gradient computed and add to the general trend ==! 168 167 CASE ( np_EXP ) ; CALL dyn_spg_exp( kt, Kmm, puu, pvv, Krhs ) ! explicit 169 CASE ( np_TS ) ; CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa , k_only_ADV) ! time-splitting168 CASE ( np_TS ) ; CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting 170 169 END SELECT 171 170 ! -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynspg_ts.F90
r14225 r14547 68 68 PUBLIC dyn_spg_ts ! called by dyn_spg 69 69 PUBLIC dyn_spg_ts_init ! - - dyn_spg_init 70 PUBLIC dyn_drg_init ! called by rk3ssh 70 71 71 72 !! Time filtered arrays at baroclinic time step: … … 80 81 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 81 82 83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshe_rhs ! RHS of ssh Eq. 84 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ue_rhs, Ve_rhs ! RHS of barotropic velocity Eq. 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: CdU_u, CdU_v ! top/bottom stress at u- & v-points 86 82 87 REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios 83 88 REAL(wp) :: r1_8 = 0.125_wp ! … … 117 122 118 123 119 SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa , k_only_ADV)124 SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 120 125 !!---------------------------------------------------------------------- 121 126 !! … … 145 150 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 146 151 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 147 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels 148 INTEGER , OPTIONAL , INTENT( in ) :: k_only_ADV ! only Advection in the RHS 152 REAL(wp), DIMENSION(jpi,jpj ,jpt), INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels 149 153 ! 150 154 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 152 156 LOGICAL :: ll_init ! =T : special startup of 2d equations 153 157 INTEGER :: noffset ! local integers : time offset for bdy update 154 REAL(wp) :: r1_Dt_b, z1_hu, z1_hv! local scalars155 REAL(wp) :: za0, za1, za2, za3 158 REAL(wp) :: z1_hu, z1_hv ! local scalars 159 REAL(wp) :: za0, za1, za2, za3 ! - - 156 160 REAL(wp) :: zztmp, zldg ! - - 157 REAL(wp) :: zhu_bck, zhv_bck, zhdiv 158 REAL(wp) :: zun_save, zvn_save 161 REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - - 162 REAL(wp) :: zun_save, zvn_save ! - - 159 163 REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 160 164 REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg … … 163 167 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 164 168 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 165 !!st#if defined key_qco166 !!st REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v167 !!st#endif168 169 ! 169 170 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 175 176 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 176 177 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 178 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace 177 179 REAL(wp) :: zt0substep ! Time of day at the beginning of the time substep 178 180 !!---------------------------------------------------------------------- … … 185 187 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 186 188 ! ! inverse of baroclinic time step 187 r1_Dt_b = 1._wp / rDt188 189 ! 189 190 ll_init = ln_bt_av ! if no time averaging, then no specific restart … … 228 229 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 229 230 ! ----------------------------------------------------------------------------- 231 #if defined key_RK3 232 ! !========================================! 233 ! !== Phase 1 for RK3 time integration ==! 234 ! !========================================! 235 ! 236 ! ! Currently, RK3 requires the forward mode and NO time filtering of barotropic variables 237 IF( kt == nit000 ) THEN 238 IF( .NOT.ln_bt_fw .OR. ln_bt_av ) CALL ctl_stop( 'dyn_spg_ts: RK3 requires ln_bt_fw=T AND ln_bt_av=F') 239 ENDIF 240 ! 241 ! ! set values computed in RK3_ssh 242 zssh_frc(:,:) = sshe_rhs(:,:) 243 zu_frc(:,:) = Ue_rhs(:,:) 244 zv_frc(:,:) = Ve_rhs(:,:) 245 zCdU_u (:,:) = CdU_u (:,:) 246 zCdU_v (:,:) = CdU_v (:,:) 247 248 !!gm ==>>> !!ts ISSUe her on en discute changer les cas ENS ENE et triad ? 249 IF( kt == nit000 .AND. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient 250 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 251 ! 252 253 #else 254 ! !========================================! 255 ! !== Phase 1 for MLF time integration ==! 256 ! !========================================! 257 ! 230 258 ! 231 259 ! 232 260 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 233 261 ! ! --------------------------- ! 234 # if defined key_qco262 # if defined key_qco 235 263 zu_frc(:,:) = SUM( e3u_0(:,:,: ) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:) 236 264 zv_frc(:,:) = SUM( e3v_0(:,:,: ) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:) 237 # else265 # else 238 266 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu(:,:,Kmm) 239 267 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv(:,:,Kmm) 240 # endif268 # endif 241 269 ! 242 270 ! … … 256 284 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 257 285 ! 258 IF( .NOT. PRESENT(k_only_ADV) ) THEN !* remove the 2D Coriolis trend 259 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 260 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 261 ! 262 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 263 & zu_trd, zv_trd ) ! ==>> out 264 ! 265 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend 266 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 267 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 268 END_2D 269 ENDIF 286 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 287 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 288 ! 289 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 290 & zu_trd, zv_trd ) ! ==>> out 291 ! 292 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend 293 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 294 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 295 END_2D 270 296 ! 271 297 ! != Add bottom stress contribution from baroclinic velocities =! 272 298 ! ! ----------------------------------------------------------- ! 273 IF( PRESENT(k_only_ADV) ) THEN !* only Advection in the RHS : provide the barotropic bottom drag coefficients 274 DO_2D( 0, 0, 0, 0 ) 275 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 276 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 277 END_2D 278 ELSE !* remove baroclinic drag AND provide the barotropic drag coefficients 279 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) 280 ENDIF 299 CALL dyn_drg_init( Kbb, Kmm, puu , pvv , puu_b , pvv_b , & ! <<= IN 300 & zu_frc, zv_frc, zCdU_u, zCdU_v ) ! =>> OUT 281 301 ! 282 302 ! != Add atmospheric pressure forcing =! … … 348 368 END IF 349 369 ! 350 # if defined key_asminc370 # if defined key_asminc 351 371 ! != Add the IAU weighted SSH increment =! 352 372 ! ! ------------------------------------ ! … … 354 374 zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 355 375 ENDIF 376 # endif 377 378 ! !== END of Phase 1 for MLF time integration ==! 356 379 #endif 380 381 357 382 ! != Fill boundary data arrays for AGRIF 358 383 ! ! ------------------------------------ … … 701 726 sshb_e (:,:) = sshn_e(:,:) 702 727 sshn_e (:,:) = ssha_e(:,:) 703 704 ! !* Sum over whole bt loop 728 ! 729 ! !* Sum over whole bt loop (except in weight average) 705 730 ! ! ---------------------- 706 za1 = wgtbtp1(jn) 707 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 708 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 709 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 710 ELSE ! Sum transports 711 IF ( .NOT.ln_wd_dl ) THEN 712 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) 713 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) 714 ELSE 715 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 716 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 717 END IF 718 ENDIF 719 ! ! Sum sea level 720 pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:) 721 731 IF( ln_bt_av ) THEN 732 za1 = wgtbtp1(jn) 733 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 734 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 735 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 736 ELSE ! Sum transports 737 IF ( .NOT.ln_wd_dl ) THEN 738 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) 739 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) 740 ELSE 741 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 742 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 743 ENDIF 744 ENDIF 745 ! ! Sum sea level 746 pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:) 747 ENDIF 748 ! 722 749 ! ! ==================== ! 723 750 END DO ! end loop ! 724 751 ! ! ==================== ! 752 753 725 754 ! ----------------------------------------------------------------------------- 726 755 ! Phase 3. update the general trend with the barotropic trend 727 756 ! ----------------------------------------------------------------------------- 728 757 ! 729 ! Set advection velocity correction: 730 IF (ln_bt_fw) THEN 758 IF(.NOT.ln_bt_av ) THEN !* Update Kaa barotropic external mode 759 puu_b(:,:,Kaa) = ua_e (:,:) 760 pvv_b(:,:,Kaa) = va_e (:,:) 761 pssh (:,:,Kaa) = ssha_e(:,:) 762 ENDIF 763 764 #if defined key_RK3 765 ! !* RK3 case 766 ! 767 ! advective transport from N to N+1 (i.e. Kbb to Kaa) 768 ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation (NOT USED) 769 vb2_b(:,:) = vn_adv(:,:) 770 ! 771 IF( iom_use("ubar") ) THEN ! RK3 single first: hu[N+1/2] = 1/2 ( hu[N] + hu[N+1] ) 772 ALLOCATE( z2d(jpi,jpj) ) 773 z2d(:,:) = 2._wp / ( hu_e(:,:) + hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) 774 CALL iom_put( "ubar", un_adv(:,:)*z2d(:,:) ) ! barotropic i-current 775 z2d(:,:) = 2._wp / ( hv_e(:,:) + hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 776 CALL iom_put( "vbar", vn_adv(:,:)*z2d(:,:) ) ! barotropic i-current 777 DEALLOCATE( z2d ) 778 ENDIF 779 ! 780 ! !== END Phase 3 for RK3 (forward mode) ==! 781 782 #else 783 ! !* MLF case 784 ! 785 ! Set advective velocity correction: 786 IF( ln_bt_fw ) THEN 731 787 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 732 788 DO_2D( 1, 1, 1, 1 ) … … 748 804 ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation 749 805 vb2_b(:,:) = vn_adv(:,:) 750 END IF 751 ENDIF 752 753 806 ENDIF 807 ENDIF 754 808 ! 755 809 ! Update barotropic trend: 756 810 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 757 811 DO jk=1,jpkm1 758 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt _b759 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt _b812 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt 813 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt 760 814 END DO 761 815 ELSE 762 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 763 #if defined key_qcoTest_FluxForm 764 ! ! 'key_qcoTest_FluxForm' : simple ssh average 765 DO_2D( 1, 0, 1, 0 ) 766 zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj) 767 zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj) 768 END_2D 769 #else 770 DO_2D( 1, 0, 1, 0 ) 771 zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 772 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 773 zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 774 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 775 END_2D 776 #endif 777 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 778 ! 779 DO jk=1,jpkm1 780 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) & 781 & * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 782 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) & 783 & * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 784 END DO 785 ! Save barotropic velocities not transport: 786 puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 787 pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 816 IF(.NOT.ln_bt_av ) THEN ! (puu_b,pvv_b)_Kaa is a velocity (hu,hv)_Kaa = (hu_e,hv_e) 817 ! 818 DO jk=1,jpkm1 819 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) & 820 & * ( puu_b(:,:,Kaa)*hu_e(:,:) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt 821 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) & 822 & * ( pvv_b(:,:,Kaa)*hv_e(:,:) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt 823 END DO 824 ! 825 ELSE ! at this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 826 ! 827 # if defined key_qcoTest_FluxForm 828 ! ! 'key_qcoTest_FluxForm' : simple ssh average 829 DO_2D( 1, 0, 1, 0 ) 830 zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj) 831 zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj) 832 END_2D 833 # else 834 DO_2D( 1, 0, 1, 0 ) 835 zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 836 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 837 zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 838 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 839 END_2D 840 # endif 841 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 842 ! 843 DO jk=1,jpkm1 844 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) & 845 & * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt 846 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) & 847 & * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt 848 END DO 849 ! Save barotropic velocities not transport: 850 puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 851 pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 852 ENDIF 788 853 ENDIF 789 854 … … 795 860 END DO 796 861 797 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN862 IF( ln_wd_dl .AND. ln_wd_dl_bc) THEN 798 863 DO jk = 1, jpkm1 799 864 puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) & 800 &+ zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk)865 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 801 866 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 802 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)867 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 803 868 END DO 804 END IF 805 806 869 ENDIF 870 807 871 CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current 808 872 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current 873 874 ! !== END Phase 3 for MLF time integration ==! 875 #endif 876 809 877 ! 810 878 #if defined key_agrif … … 834 902 END SUBROUTINE dyn_spg_ts 835 903 836 837 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)904 905 SUBROUTINE ts_wgt( ll_av, ll_fw, Kpit, zwgt1, zwgt2) 838 906 !!--------------------------------------------------------------------- 839 907 !! *** ROUTINE ts_wgt *** … … 841 909 !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 842 910 !!---------------------------------------------------------------------- 843 LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true. 844 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 845 INTEGER, INTENT(inout) :: jpit ! cycle length 846 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights 847 zwgt2 ! Secondary weights 848 849 INTEGER :: jic, jn, ji ! temporary integers 850 REAL(wp) :: za1, za2 851 !!---------------------------------------------------------------------- 852 911 LOGICAL, INTENT(in ) :: ll_av ! temporal averaging=.true. 912 LOGICAL, INTENT(in ) :: ll_fw ! forward time splitting =.true. 913 INTEGER, INTENT(inout) :: Kpit ! cycle length 914 !! 915 INTEGER :: jic, jn, ji ! local integers 916 REAL(wp) :: za1, za2 ! loca scalars 917 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, zwgt2 ! Primary & Secondary weights 918 !!---------------------------------------------------------------------- 919 ! 853 920 zwgt1(:) = 0._wp 854 921 zwgt2(:) = 0._wp 855 856 ! Set time index when averaged value is requested 857 IF (ll_fw) THEN 858 jic = nn_e 859 ELSE 860 jic = 2 * nn_e 861 ENDIF 862 863 ! Set primary weights: 864 IF (ll_av) THEN 865 ! Define simple boxcar window for primary weights 866 ! (width = nn_e, centered around jic) 867 SELECT CASE ( nn_bt_flt ) 868 CASE( 0 ) ! No averaging 869 zwgt1(jic) = 1._wp 870 jpit = jic 871 872 CASE( 1 ) ! Boxcar, width = nn_e 873 DO jn = 1, 3*nn_e 874 za1 = ABS(float(jn-jic))/float(nn_e) 875 IF (za1 < 0.5_wp) THEN 876 zwgt1(jn) = 1._wp 877 jpit = jn 878 ENDIF 879 ENDDO 880 881 CASE( 2 ) ! Boxcar, width = 2 * nn_e 882 DO jn = 1, 3*nn_e 883 za1 = ABS(float(jn-jic))/float(nn_e) 884 IF (za1 < 1._wp) THEN 885 zwgt1(jn) = 1._wp 886 jpit = jn 887 ENDIF 888 ENDDO 889 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 922 ! 923 ! !== Set time index when averaged value is requested ==! 924 IF (ll_fw) THEN ; jic = nn_e 925 ELSE ; jic = 2 * nn_e 926 ENDIF 927 ! 928 ! !== Set primary weights ==! 929 ! 930 IF (ll_av) THEN != Define simple boxcar window for primary weights 931 ! ! (width = nn_e, centered around jic) 932 SELECT CASE( nn_bt_flt ) 933 ! 934 CASE( 0 ) ! No averaging 935 zwgt1(jic) = 1._wp 936 Kpit = jic 937 ! 938 CASE( 1 ) ! Boxcar, width = nn_e 939 DO jn = 1, 3*nn_e 940 za1 = ABS( REAL( jn-jic, wp) ) / REAL( nn_e, wp ) 941 IF( za1 < 0.5_wp ) THEN 942 zwgt1(jn) = 1._wp 943 Kpit = jn 944 ENDIF 945 END DO 946 ! 947 CASE( 2 ) ! Boxcar, width = 2 * nn_e 948 DO jn = 1, 3*nn_e 949 za1 = ABS(REAL( jn-jic, wp) ) / REAL( nn_e, wp ) 950 IF( za1 < 1._wp ) THEN 951 zwgt1(jn) = 1._wp 952 Kpit = jn 953 ENDIF 954 END DO 955 ! 956 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 957 ! 890 958 END SELECT 891 959 892 ELSE !No time averaging960 ELSE != No time averaging 893 961 zwgt1(jic) = 1._wp 894 jpit = jic895 ENDIF 896 897 ! Set secondary weights898 DO jn = 1, jpit899 DO ji = jn, jpit900 901 END DO962 Kpit = jic 963 ENDIF 964 ! 965 ! !== Set secondary weights ==! 966 DO jn = 1, Kpit 967 DO ji = jn, Kpit 968 zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 969 END DO 902 970 END DO 903 904 ! Normalize weigths:905 za1 = 1._wp / SUM( zwgt1(1:jpit))906 za2 = 1._wp / SUM( zwgt2(1:jpit))907 DO jn = 1, jpit908 zwgt1(jn) = zwgt1(jn) * za1909 zwgt2(jn) = zwgt2(jn) * za2971 ! 972 ! !== Normalize weigths ==! 973 za1 = 1._wp / SUM( zwgt1(1:Kpit) ) 974 za2 = 1._wp / SUM( zwgt2(1:Kpit) ) 975 DO jn = 1, Kpit 976 zwgt1(jn) = zwgt1(jn) * za1 977 zwgt2(jn) = zwgt2(jn) * za2 910 978 END DO 911 979 ! … … 1023 1091 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e 1024 1092 ENDIF 1025 1093 ! 1026 1094 IF(ln_bt_av) THEN 1027 1095 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on ' … … 1029 1097 IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables ' 1030 1098 ENDIF 1031 !1032 1099 ! 1033 1100 IF(ln_bt_fw) THEN … … 1088 1155 !! and update depths at T- points (ht) at each barotropic time step 1089 1156 !! 1090 !! Compute zwz = f / ( height of the water colomn ) 1157 !! Compute zwz = f/h (potential planetary voricity) 1158 !! Compute ftne, ftnw, ftse, ftsw (triad of potential planetary voricity) 1091 1159 !!---------------------------------------------------------------------- 1092 1160 INTEGER, INTENT(in) :: Kmm ! Time index … … 1138 1206 ! 1139 1207 END SELECT 1140 1208 ! 1141 1209 END SUBROUTINE dyn_cor_2d_init 1142 1210 … … 1287 1355 !! ** Purpose : 1288 1356 !!---------------------------------------------------------------------- 1357 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn 1358 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1359 ! 1289 1360 INTEGER :: ji ,jj ! dummy loop indices 1290 1361 LOGICAL :: ll_tmp1, ll_tmp2 1291 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn 1292 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1293 !!---------------------------------------------------------------------- 1362 !!---------------------------------------------------------------------- 1363 ! 1294 1364 DO_2D( 0, 0, 0, 0 ) 1295 1365 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > & … … 1330 1400 ENDIF 1331 1401 END_2D 1332 1402 ! 1333 1403 END SUBROUTINE wad_spg 1334 1404 … … 1374 1444 ! 1375 1445 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1376 1377 1446 DO_2D( 0, 0, 0, 0 ) 1378 1447 ikbu = mbku(ji,jj) … … 1381 1450 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 1382 1451 END_2D 1383 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1384 1452 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1385 1453 DO_2D( 0, 0, 0, 0 ) 1386 1454 ikbu = mbku(ji,jj) … … 1400 1468 END_2D 1401 1469 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1402 1403 1470 DO_2D( 0, 0, 0, 0 ) 1404 1471 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) … … 1412 1479 ! 1413 1480 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1414 1415 1481 DO_2D( 0, 0, 0, 0 ) 1416 1482 iktu = miku(ji,jj) … … 1419 1485 zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 1420 1486 END_2D 1421 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1422 1487 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1423 1488 DO_2D( 0, 0, 0, 0 ) 1424 1489 iktu = miku(ji,jj) … … 1428 1493 END_2D 1429 1494 ENDIF 1430 ! 1431 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1432 1495 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1433 1496 DO_2D( 0, 0, 0, 0 ) 1434 1497 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) … … 1439 1502 ! 1440 1503 END SUBROUTINE dyn_drg_init 1504 1441 1505 1442 1506 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in … … 1477 1541 END SUBROUTINE ts_bck_interp 1478 1542 1479 1480 1543 !!====================================================================== 1481 1544 END MODULE dynspg_ts -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90
r14233 r14547 51 51 IMPLICIT NONE 52 52 PRIVATE 53 54 INTERFACE dyn_vor 55 MODULE PROCEDURE dyn_vor_3D, dyn_vor_RK3 56 END INTERFACE 53 57 54 58 PUBLIC dyn_vor ! routine called by step.F90 … … 74 78 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 75 79 76 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 77 ! ! associated indices: 80 ! !: choice of calculated vorticity 81 INTEGER, PUBLIC :: ncor, nrvm, ntot ! Coriolis, relative vorticity, total vorticity 82 ! ! associated indices: 78 83 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) 79 84 INTEGER, PUBLIC, PARAMETER :: np_RVO = 2 ! relative vorticity … … 103 108 !!---------------------------------------------------------------------- 104 109 CONTAINS 105 106 SUBROUTINE dyn_vor ( kt, Kmm, puu, pvv, Krhs )110 111 SUBROUTINE dyn_vor_3D( kt, Kmm, puu, pvv, Krhs ) 107 112 !!---------------------------------------------------------------------- 108 113 !! … … 121 126 !!---------------------------------------------------------------------- 122 127 ! 123 IF( ln_timing ) CALL timing_start('dyn_vor ')128 IF( ln_timing ) CALL timing_start('dyn_vor_3D') 124 129 ! 125 130 IF( l_trddyn ) THEN !== trend diagnostics case : split the added trend in two parts ==! … … 209 214 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 210 215 ! 211 IF( ln_timing ) CALL timing_stop('dyn_vor') 212 ! 213 END SUBROUTINE dyn_vor 216 IF( ln_timing ) CALL timing_stop('dyn_vor_3D') 217 ! 218 END SUBROUTINE dyn_vor_3D 219 220 221 SUBROUTINE dyn_vor_RK3( kt, Kmm, puu, pvv, Krhs, knoco ) 222 !!---------------------------------------------------------------------- 223 !! 224 !! ** Purpose : compute the lateral ocean tracer physics. 225 !! 226 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 227 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative 228 !! and planetary vorticity trends) and send them to trd_dyn 229 !! for futher diagnostics (l_trddyn=T) 230 !!---------------------------------------------------------------------- 231 INTEGER , INTENT(in ) :: kt ! ocean time-step index 232 INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices 233 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation 234 INTEGER , INTENT(in ) :: knoco ! specified vorticity trend (= np_MET or np_RVO) 235 !!---------------------------------------------------------------------- 236 ! 237 IF( ln_timing ) CALL timing_start('dyn_vor_RK3') 238 ! 239 ! !== total vorticity trend added to the general trend ==! 240 !!st WARNING 22/02 !!!!!!!! stoke drift or not stoke drift ? => bar to do later !!! 241 !! stoke drift a garder pas vortex force a priori !! 242 !! ATTENTION déja appelé avec Kbb !! 243 244 ! 245 SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==! 246 CASE( np_ENT ) !* energy conserving scheme (T-pts) 247 CALL vor_enT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 248 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 249 CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 250 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 251 CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 252 ENDIF 253 CASE( np_EET ) !* energy conserving scheme (een scheme using e3t) 254 CALL vor_eeT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 255 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 256 CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 257 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 258 CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 259 ENDIF 260 CASE( np_ENE ) !* energy conserving scheme 261 CALL vor_ene( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 262 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 263 CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 264 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 265 CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 266 ENDIF 267 CASE( np_ENS ) !* enstrophy conserving scheme 268 CALL vor_ens( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 269 270 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 271 CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 272 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 273 CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 274 ENDIF 275 CASE( np_MIX ) !* mixed ene-ens scheme 276 CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! relative vorticity or metric trend (ens) 277 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 278 IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add vortex force 279 CASE( np_EEN ) !* energy and enstrophy conserving scheme 280 CALL vor_een( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 281 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 282 CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 283 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 284 CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 285 ENDIF 286 END SELECT 287 ! 288 ! ! print sum trends (used for debugging) 289 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, & 290 & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 291 ! 292 IF( ln_timing ) CALL timing_stop('dyn_vor_RK3') 293 ! 294 END SUBROUTINE dyn_vor_RK3 214 295 215 296 -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynzdf.F90
r13497 r14547 70 70 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 71 71 ! 72 INTEGER :: ji, jj, jk ! dummy loop indices73 INTEGER :: iku, ikv ! local integers74 REAL(wp) :: zzwi, ze3ua, z dt! local scalars75 REAL(wp) :: zzws, ze3va ! - -76 REAL(wp) :: z1_e3ua, z1_e3va ! - -77 REAL(wp) :: zWu , zWv ! - -78 REAL(wp) :: zWui, zWvi ! - -79 REAL(wp) :: zWus, zWvs ! - -72 INTEGER :: ji, jj, jk ! dummy loop indices 73 INTEGER :: iku, ikv ! local integers 74 REAL(wp) :: zzwi, ze3ua, zDt_2 ! local scalars 75 REAL(wp) :: zzws, ze3va ! - - 76 REAL(wp) :: z1_e3ua, z1_e3va ! - - 77 REAL(wp) :: zWu , zWv ! - - 78 REAL(wp) :: zWui, zWvi ! - - 79 REAL(wp) :: zWus, zWvs ! - - 80 80 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace 81 81 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - … … 93 93 ENDIF 94 94 ENDIF 95 ! 96 zDt_2 = rDt * 0.5_wp 97 ! 95 98 ! !* explicit top/bottom drag case 96 99 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa)) … … 138 141 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 139 142 & + r_vvl * e3v(ji,jj,ikv,Kaa) 140 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua141 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va143 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 144 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 142 145 END_2D 143 146 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF) … … 149 152 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 150 153 & + r_vvl * e3v(ji,jj,ikv,Kaa) 151 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua152 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va154 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 155 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 153 156 END_2D 154 157 END IF … … 158 161 ! 159 162 ! !* Matrix construction 160 zdt = rDt * 0.5161 163 IF( ln_zad_Aimp ) THEN !! 162 164 SELECT CASE( nldf_dyn ) … … 165 167 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 166 168 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 167 zzwi = - z dt* ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &169 zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 168 170 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 169 zzws = - z dt* ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &171 zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 170 172 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 171 173 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 172 174 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 173 zwi(ji,jj,jk) = zzwi + z dt* MIN( zWui, 0._wp )174 zws(ji,jj,jk) = zzws - z dt* MAX( zWus, 0._wp )175 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + z dt* ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )175 zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp ) 176 zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp ) 177 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 176 178 END_3D 177 179 CASE DEFAULT ! iso-level lateral mixing … … 179 181 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & ! after scale factor at U-point 180 182 & + r_vvl * e3u(ji,jj,jk,Kaa) 181 zzwi = - z dt* ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &183 zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 182 184 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 183 zzws = - z dt* ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &185 zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 184 186 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 185 187 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 186 188 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 187 zwi(ji,jj,jk) = zzwi + z dt* MIN( zWui, 0._wp )188 zws(ji,jj,jk) = zzws - z dt* MAX( zWus, 0._wp )189 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + z dt* ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )189 zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp ) 190 zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp ) 191 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 190 192 END_3D 191 193 END SELECT … … 194 196 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 195 197 & + r_vvl * e3u(ji,jj,1,Kaa) 196 zzws = - z dt* ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) &198 zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) & 197 199 & / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 198 200 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 199 zws(ji,jj,1 ) = zzws - z dt* MAX( zWus, 0._wp )200 zwd(ji,jj,1 ) = 1._wp - zzws - z dt* ( MIN( zWus, 0._wp ) )201 zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWus, 0._wp ) 202 zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) ) 201 203 END_2D 202 204 ELSE … … 206 208 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 207 209 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 208 zzwi = - z dt* ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &210 zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 209 211 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 210 zzws = - z dt* ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &212 zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 211 213 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 212 214 zwi(ji,jj,jk) = zzwi … … 218 220 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 219 221 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 220 zzwi = - z dt* ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &222 zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 221 223 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 222 zzws = - z dt* ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &224 zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 223 225 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 224 226 zwi(ji,jj,jk) = zzwi … … 245 247 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 246 248 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 247 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua249 zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 248 250 END_2D 249 251 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit) … … 253 255 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 254 256 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 255 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua257 zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 256 258 END_2D 257 259 END IF … … 280 282 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 281 283 & + r_vvl * e3u(ji,jj,1,Kaa) 282 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp* ( utau_b(ji,jj) + utau(ji,jj) ) &284 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utau(ji,jj) ) & 283 285 & / ( ze3ua * rho0 ) * umask(ji,jj,1) 284 286 END_2D … … 297 299 ! 298 300 ! !* Matrix construction 299 zdt = rDt * 0.5300 301 IF( ln_zad_Aimp ) THEN !! 301 302 SELECT CASE( nldf_dyn ) … … 304 305 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 305 306 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 306 zzwi = - z dt* ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &307 zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 307 308 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 308 zzws = - z dt* ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &309 zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 309 310 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 310 311 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 311 312 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 312 zwi(ji,jj,jk) = zzwi + z dt* MIN( zWvi, 0._wp )313 zws(ji,jj,jk) = zzws - z dt* MAX( zWvs, 0._wp )314 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - z dt* ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )313 zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp ) 314 zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp ) 315 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 315 316 END_3D 316 317 CASE DEFAULT ! iso-level lateral mixing … … 318 319 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 319 320 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 zzwi = - z dt* ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &321 zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 321 322 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 322 zzws = - z dt* ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &323 zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 323 324 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 324 325 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 325 326 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 326 zwi(ji,jj,jk) = zzwi + z dt* MIN( zWvi, 0._wp )327 zws(ji,jj,jk) = zzws - z dt* MAX( zWvs, 0._wp )328 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - z dt* ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )327 zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp ) 328 zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp ) 329 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 329 330 END_3D 330 331 END SELECT … … 333 334 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 334 335 & + r_vvl * e3v(ji,jj,1,Kaa) 335 zzws = - z dt* ( avm(ji,jj+1,2) + avm(ji,jj,2) ) &336 zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) & 336 337 & / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 337 338 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 338 zws(ji,jj,1 ) = zzws - z dt* MAX( zWvs, 0._wp )339 zwd(ji,jj,1 ) = 1._wp - zzws - z dt* ( MIN( zWvs, 0._wp ) )339 zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp ) 340 zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) ) 340 341 END_2D 341 342 ELSE … … 345 346 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 346 347 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 347 zzwi = - z dt* ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &348 zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 348 349 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 349 zzws = - z dt* ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &350 zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 350 351 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 351 352 zwi(ji,jj,jk) = zzwi … … 357 358 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 358 359 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 359 zzwi = - z dt* ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &360 zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 360 361 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 361 zzws = - z dt* ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &362 zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 362 363 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 363 364 zwi(ji,jj,jk) = zzwi … … 383 384 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 384 385 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 385 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va386 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 386 387 END_2D 387 388 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN … … 390 391 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 391 392 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 392 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va393 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 393 394 END_2D 394 395 ENDIF … … 417 418 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 418 419 & + r_vvl * e3v(ji,jj,1,Kaa) 419 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp *( vtau_b(ji,jj) + vtau(ji,jj) ) &420 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtau(ji,jj) ) & 420 421 & / ( ze3va * rho0 ) * vmask(ji,jj,1) 421 422 END_2D … … 432 433 ! 433 434 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 434 ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:)435 ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:)435 ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:) 436 ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) )*r1_Dt - ztrdv(:,:,:) 436 437 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 437 438 DEALLOCATE( ztrdu, ztrdv ) -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/sshwzv.F90
r14205 r14547 17 17 !! ssh_nxt : after ssh 18 18 !! ssh_atf : time filter the ssh arrays 19 !! wzv : compute now vertical velocity 19 !! wzv : generic interface of vertical velocity calculation 20 !! wzv_MLF : MLF: compute NOW vertical velocity 21 !! wzv_RK3 : RK3: Compute a vertical velocity 20 22 !!---------------------------------------------------------------------- 21 23 USE oce ! ocean dynamics and tracers variables … … 44 46 IMPLICIT NONE 45 47 PRIVATE 48 ! !! * Interface 49 INTERFACE wzv 50 MODULE PROCEDURE wzv_MLF, wzv_RK3 51 END INTERFACE 46 52 47 53 PUBLIC ssh_nxt ! called by step.F90 … … 77 83 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index 78 84 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height 79 ! 85 ! 80 86 INTEGER :: jk ! dummy loop index 81 87 REAL(wp) :: zcoef ! local scalar 82 88 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 89 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace 83 90 !!---------------------------------------------------------------------- 84 91 ! … … 92 99 ! 93 100 zcoef = 0.5_wp * r1_rho0 94 95 101 ! !------------------------------! 96 102 ! ! After Sea Surface Height ! 97 103 ! !------------------------------! 98 104 IF(ln_wd_il) THEN 99 CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )105 CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv ) 100 106 ENDIF 101 107 … … 133 139 END SUBROUTINE ssh_nxt 134 140 135 136 SUBROUTINE wzv ( kt, Kbb, Kmm, Kaa, pww )137 !!---------------------------------------------------------------------- 138 !! *** ROUTINE wzv ***141 142 SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww ) 143 !!---------------------------------------------------------------------- 144 !! *** ROUTINE wzv_MLF *** 139 145 !! 140 146 !! ** Purpose : compute the now vertical velocity … … 157 163 !!---------------------------------------------------------------------- 158 164 ! 159 IF( ln_timing ) CALL timing_start('wzv ')165 IF( ln_timing ) CALL timing_start('wzv_MLF') 160 166 ! 161 167 IF( kt == nit000 ) THEN 162 168 IF(lwp) WRITE(numout,*) 163 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '164 IF(lwp) WRITE(numout,*) '~~~~~ 169 IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity ' 170 IF(lwp) WRITE(numout,*) '~~~~~~~' 165 171 ! 166 172 pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) … … 193 199 END DO 194 200 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 195 DEALLOCATE( zhdiv ) 201 DEALLOCATE( zhdiv ) 196 202 ! !=================================! 197 203 ELSEIF( ln_linssh ) THEN !== linear free surface cases ==! … … 204 210 ! !==========================================! 205 211 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 206 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 207 & + r1_Dt * ( e3t(:,:,jk,Kaa) & 208 & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 212 ! ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]] 213 !!gm old 214 ! pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 215 ! & + r1_Dt * ( e3t(:,:,jk,Kaa) & 216 ! & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 217 !!gm slightly faster : 218 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 219 & + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) ) ) * tmask(:,:,jk) 220 !!gm end 221 222 209 223 END DO 210 224 ENDIF … … 256 270 #endif 257 271 ! 258 IF( ln_timing ) CALL timing_stop('wzv') 259 ! 260 END SUBROUTINE wzv 272 IF( ln_timing ) CALL timing_stop('wzv_MLF') 273 ! 274 END SUBROUTINE wzv_MLF 275 276 277 SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww ) 278 !!---------------------------------------------------------------------- 279 !! *** ROUTINE wzv_RK3 *** 280 !! 281 !! ** Purpose : compute the now vertical velocity 282 !! 283 !! ** Method : - Using the incompressibility hypothesis, the vertical 284 !! velocity is computed by integrating the horizontal divergence 285 !! from the bottom to the surface minus the scale factor evolution. 286 !! The boundary conditions are w=0 at the bottom (no flux) and. 287 !! 288 !! ** action : pww : now vertical velocity 289 !! 290 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 291 !!---------------------------------------------------------------------- 292 INTEGER , INTENT(in ) :: kt ! time step 293 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices 294 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: puu, pvv ! horizontal velocity at Kmm 295 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm 296 ! 297 INTEGER :: ji, jj, jk ! dummy loop indices 298 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv 299 REAL(wp) , DIMENSION(jpi,jpj,jpk) :: ze3div 300 !!---------------------------------------------------------------------- 301 ! 302 IF( ln_timing ) CALL timing_start('wzv_RK3') 303 ! 304 IF( kt == nit000 ) THEN 305 IF(lwp) WRITE(numout,*) 306 IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity ' 307 IF(lwp) WRITE(numout,*) '~~~~~ ' 308 ! 309 pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 310 ENDIF 311 ! 312 CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div ) 313 ! !------------------------------! 314 ! ! Now Vertical Velocity ! 315 ! !------------------------------! 316 ! 317 ! !===============================! 318 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==! 319 ! !===============================! 320 ALLOCATE( zhdiv(jpi,jpj,jpk) ) 321 ! 322 DO jk = 1, jpkm1 323 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 324 ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 325 DO_2D( 0, 0, 0, 0 ) 326 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) ) 327 END_2D 328 END DO 329 CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 330 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 331 ! ! Same question holds for hdiv. Perhaps just for security 332 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 333 ! computation of w 334 pww(:,:,jk) = pww(:,:,jk+1) - ( ze3div(:,:,jk) + zhdiv(:,:,jk) & 335 & + r1_Dt * ( e3t(:,:,jk,Kaa) & 336 & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 337 END DO 338 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 339 DEALLOCATE( zhdiv ) 340 ! !=================================! 341 ELSEIF( ln_linssh ) THEN !== linear free surface cases ==! 342 ! !=================================! 343 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 344 pww(:,:,jk) = pww(:,:,jk+1) - ze3div(:,:,jk) 345 END DO 346 ! !==========================================! 347 ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco') 348 ! !==========================================! 349 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 350 ! ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]] 351 pww(:,:,jk) = pww(:,:,jk+1) - ( ze3div(:,:,jk) & 352 & + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) ) ) * tmask(:,:,jk) 353 END DO 354 ENDIF 355 356 IF( ln_bdy ) THEN 357 DO jk = 1, jpkm1 358 pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:) 359 END DO 360 ENDIF 361 ! 362 #if defined key_agrif 363 IF( .NOT. AGRIF_Root() ) THEN 364 ! 365 ! Mask vertical velocity at first/last columns/row 366 ! inside computational domain (cosmetic) 367 DO jk = 1, jpkm1 368 IF( lk_west ) THEN ! --- West --- ! 369 DO ji = mi0(2+nn_hls), mi1(2+nn_hls) 370 DO jj = 1, jpj 371 pww(ji,jj,jk) = 0._wp 372 END DO 373 END DO 374 ENDIF 375 IF( lk_east ) THEN ! --- East --- ! 376 DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls) 377 DO jj = 1, jpj 378 pww(ji,jj,jk) = 0._wp 379 END DO 380 END DO 381 ENDIF 382 IF( lk_south ) THEN ! --- South --- ! 383 DO jj = mj0(2+nn_hls), mj1(2+nn_hls) 384 DO ji = 1, jpi 385 pww(ji,jj,jk) = 0._wp 386 END DO 387 END DO 388 ENDIF 389 IF( lk_north ) THEN ! --- North --- ! 390 DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls) 391 DO ji = 1, jpi 392 pww(ji,jj,jk) = 0._wp 393 END DO 394 END DO 395 ENDIF 396 ! 397 END DO 398 ! 399 ENDIF 400 #endif 401 ! 402 IF( ln_timing ) CALL timing_stop('wzv_RK3') 403 ! 404 END SUBROUTINE wzv_RK3 261 405 262 406 -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/eosbn2.F90
r14131 r14547 56 56 ! !! * Interface 57 57 INTERFACE eos 58 MODULE PROCEDURE eos_insitu , eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d58 MODULE PROCEDURE eos_insitu_New, eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d 59 59 END INTERFACE 60 60 ! … … 188 188 !!---------------------------------------------------------------------- 189 189 CONTAINS 190 191 SUBROUTINE eos_insitu_New( pts, Knn, prd ) 192 !!---------------------------------------------------------------------- 193 !! *** ROUTINE eos_insitu *** 194 !! 195 !! ** Purpose : Compute the in situ density (ratio rho/rho0) from 196 !! potential temperature and salinity using an equation of state 197 !! selected in the nameos namelist 198 !! 199 !! ** Method : prd(t,s,z) = ( rho(t,s,z) - rho0 ) / rho0 200 !! with prd in situ density anomaly no units 201 !! t TEOS10: CT or EOS80: PT Celsius 202 !! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu 203 !! z depth meters 204 !! rho in situ density kg/m^3 205 !! rho0 reference density kg/m^3 206 !! 207 !! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z). 208 !! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg 209 !! 210 !! ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z). 211 !! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu 212 !! 213 !! ln_seos : simplified equation of state 214 !! prd(t,s,z) = ( -a0*(1+lambda/2*(T-T0)+mu*z+nu*(S-S0))*(T-T0) + b0*(S-S0) ) / rho0 215 !! linear case function of T only: rn_alpha<>0, other coefficients = 0 216 !! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0 217 !! Vallis like equation: use default values of coefficients 218 !! 219 !! ** Action : compute prd , the in situ density (no units) 220 !! 221 !! References : Roquet et al, Ocean Modelling, in preparation (2014) 222 !! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006 223 !! TEOS-10 Manual, 2010 224 !!---------------------------------------------------------------------- 225 REAL(wp), DIMENSION(:,:,:,:,:), INTENT(in ) :: pts ! T-S 226 INTEGER , INTENT(in ) :: Knn ! time-level 227 REAL(wp), DIMENSION(:,:,: ), INTENT( out) :: prd ! in situ density 228 ! 229 INTEGER :: ji, jj, jk ! dummy loop indices 230 REAL(wp) :: zt , zh , zs , ztm ! local scalars 231 REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - - 232 !!---------------------------------------------------------------------- 233 ! 234 IF( ln_timing ) CALL timing_start('eos-insitu') 235 ! 236 SELECT CASE( neos ) 237 ! 238 CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==! 239 ! 240 DO_3D(nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 241 ! 242 zh = gdept(ji,jj,jk,Knn) * r1_Z0 ! depth 243 zt = pts (ji,jj,jk,jp_tem,Knn) * r1_T0 ! temperature 244 zs = SQRT( ABS( pts(ji,jj,jk,jp_sal,Knn) + rdeltaS ) * r1_S0 ) ! square root salinity 245 ztm = tmask(ji,jj,jk) ! tmask 246 ! 247 zn3 = EOS013*zt & 248 & + EOS103*zs+EOS003 249 ! 250 zn2 = (EOS022*zt & 251 & + EOS112*zs+EOS012)*zt & 252 & + (EOS202*zs+EOS102)*zs+EOS002 253 ! 254 zn1 = (((EOS041*zt & 255 & + EOS131*zs+EOS031)*zt & 256 & + (EOS221*zs+EOS121)*zs+EOS021)*zt & 257 & + ((EOS311*zs+EOS211)*zs+EOS111)*zs+EOS011)*zt & 258 & + (((EOS401*zs+EOS301)*zs+EOS201)*zs+EOS101)*zs+EOS001 259 ! 260 zn0 = (((((EOS060*zt & 261 & + EOS150*zs+EOS050)*zt & 262 & + (EOS240*zs+EOS140)*zs+EOS040)*zt & 263 & + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt & 264 & + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt & 265 & + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt & 266 & + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000 267 ! 268 zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 269 ! 270 prd(ji,jj,jk) = ( zn * r1_rho0 - 1._wp ) * ztm ! density anomaly (masked) 271 ! 272 END_3D 273 ! 274 CASE( np_seos ) !== simplified EOS ==! 275 ! 276 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) 277 zt = pts (ji,jj,jk,jp_tem,Knn) - 10._wp 278 zs = pts (ji,jj,jk,jp_sal,Knn) - 35._wp 279 zh = gdept(ji,jj,jk,Knn) 280 ztm = tmask(ji,jj,jk) 281 ! 282 zn = - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt + rn_mu1*zh ) * zt & 283 & + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs - rn_mu2*zh ) * zs & 284 & - rn_nu * zt * zs 285 ! 286 prd(ji,jj,jk) = zn * r1_rho0 * ztm ! density anomaly (masked) 287 END_3D 288 ! 289 END SELECT 290 ! 291 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', kdim=jpk ) 292 ! 293 IF( ln_timing ) CALL timing_stop('eos-insitu') 294 ! 295 END SUBROUTINE eos_insitu_New 296 190 297 191 298 SUBROUTINE eos_insitu( pts, prd, pdep ) -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/TRA/tranpc.F90
r14215 r14547 71 71 LOGICAL :: l_bottom_reached, l_column_treated 72 72 REAL(wp) :: zta, zalfa, zsum_temp, zsum_alfa, zaw, zdz, zsum_z 73 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw , z1_rDt73 REAL(wp) :: zsa, zbeta, zsum_sali, zsum_beta, zbw, zrw 74 74 REAL(wp), PARAMETER :: zn2_zero = 1.e-14_wp ! acceptance criteria for neutrality (N2==0) 75 75 REAL(wp), DIMENSION( jpk ) :: zvn2 ! vertical profile of N2 at 1 given point... … … 311 311 ! 312 312 IF( l_trdtra ) THEN ! send the Non penetrative mixing trends for diagnostic 313 z1_rDt = 1._wp / (2._wp * rn_Dt) 314 ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * z1_rDt 315 ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * z1_rDt 313 ztrdt(:,:,:) = ( pts(:,:,:,jp_tem,Kaa) - ztrdt(:,:,:) ) * r1_Dt 314 ztrds(:,:,:) = ( pts(:,:,:,jp_sal,Kaa) - ztrds(:,:,:) ) * r1_Dt 316 315 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_npc, ztrdt ) 317 316 CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_npc, ztrds ) -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/nemogcm.F90
r14239 r14547 67 67 #endif 68 68 #if defined key_qco || defined key_linssh 69 # if defined key_RK3 70 USE stprk3 71 # else 69 72 USE stpmlf ! NEMO time-stepping (stp_MLF routine) 73 # endif 70 74 #else 71 75 USE step ! NEMO time-stepping (stp routine) … … 163 167 ! 164 168 # if defined key_qco || defined key_linssh 169 # if defined key_RK3 170 CALL stp_RK3 171 # else 165 172 CALL stp_MLF 173 # endif 166 174 # else 167 175 CALL stp … … 184 192 ! 185 193 # if defined key_qco || defined key_linssh 194 # if defined key_RK3 195 CALL stp_RK3( istp ) 196 # else 186 197 CALL stp_MLF( istp ) 198 # endif 187 199 # else 188 200 CALL stp ( istp ) -
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stpmlf.F90
r14381 r14547 62 62 # include "do_loop_substitute.h90" 63 63 # include "domzgr_substitute.h90" 64 # include "do_loop_substitute.h90"65 64 !!---------------------------------------------------------------------- 66 65 !! NEMO/OCE 4.0 , NEMO Consortium (2018)
Note: See TracChangeset
for help on using the changeset viewer.
