- Timestamp:
- 2020-06-24T14:38:26+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN
- Files:
-
- 1 added
- 13 edited
-
divhor.F90 (modified) (6 diffs)
-
dynadv_cen2.F90 (modified) (3 diffs)
-
dynadv_ubs.F90 (modified) (3 diffs)
-
dynatf.F90 (modified) (15 diffs)
-
dynatf_qco.F90 (added)
-
dynhpg.F90 (modified) (26 diffs)
-
dynldf_iso.F90 (modified) (14 diffs)
-
dynldf_lap_blp.F90 (modified) (5 diffs)
-
dynspg_ts.F90 (modified) (10 diffs)
-
dynvor.F90 (modified) (22 diffs)
-
dynzad.F90 (modified) (2 diffs)
-
dynzdf.F90 (modified) (21 diffs)
-
sshwzv.F90 (modified) (23 diffs)
-
wet_dry.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/divhor.F90
r12377 r13151 21 21 USE dom_oce ! ocean space and time domain 22 22 USE sbc_oce, ONLY : ln_rnf ! river runoff 23 USE sbcrnf , ONLY : sbc_rnf_div ! river runoff 23 USE sbcrnf , ONLY : sbc_rnf_div ! river runoff 24 24 USE isf_oce, ONLY : ln_isf ! ice shelf 25 25 USE isfhdiv, ONLY : isf_hdiv ! ice shelf 26 #if defined key_asminc 26 #if defined key_asminc 27 27 USE asminc ! Assimilation increment 28 28 #endif … … 40 40 !! * Substitutions 41 41 # include "do_loop_substitute.h90" 42 # include "domzgr_substitute.h90" 42 43 !!---------------------------------------------------------------------- 43 44 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 44 !! $Id$ 45 !! $Id$ 45 46 !! Software governed by the CeCILL license (see ./LICENSE) 46 47 !!---------------------------------------------------------------------- … … 50 51 !!---------------------------------------------------------------------- 51 52 !! *** ROUTINE div_hor *** 52 !! 53 !! 53 54 !! ** Purpose : compute the horizontal divergence at now time-step 54 55 !! 55 56 !! ** Method : the now divergence is computed as : 56 57 !! hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 57 !! and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 58 !! and correct with runoff inflow (div_rnf) and cross land flow (div_cla) 58 59 !! 59 60 !! ** Action : - update hdiv, the now horizontal divergence … … 78 79 DO_3D_00_00( 1, jpkm1 ) 79 80 hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * uu(ji ,jj,jk,Kmm) & 80 & 81 & 82 & 81 & - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * uu(ji-1,jj,jk,Kmm) & 82 & + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vv(ji,jj ,jk,Kmm) & 83 & - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vv(ji,jj-1,jk,Kmm) ) & 83 84 & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) 84 85 END_3D … … 95 96 IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== runoffs ==! (update hdiv field) 96 97 ! 97 #if defined key_asminc 98 #if defined key_asminc 98 99 IF( ln_sshinc .AND. ln_asmiau ) CALL ssh_asm_div( kt, Kbb, Kmm, hdiv ) !== SSH assimilation ==! (update hdiv field) 99 ! 100 ! 100 101 #endif 101 102 ! … … 107 108 ! 108 109 END SUBROUTINE div_hor 109 110 110 111 !!====================================================================== 111 112 END MODULE divhor -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynadv_cen2.F90
r12377 r13151 28 28 !! * Substitutions 29 29 # include "do_loop_substitute.h90" 30 # include "domzgr_substitute.h90" 30 31 !!---------------------------------------------------------------------- 31 32 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 79 80 DO_2D_00_00 80 81 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 81 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 82 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) & 83 & / e3u(ji,jj,jk,Kmm) 82 84 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 83 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 85 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) & 86 & / e3v(ji,jj,jk,Kmm) 84 87 END_2D 85 88 END DO … … 115 118 END DO 116 119 DO_3D_00_00( 1, jpkm1 ) 117 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 118 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 120 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) & 121 & / e3u(ji,jj,jk,Kmm) 122 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) & 123 & / e3v(ji,jj,jk,Kmm) 119 124 END_3D 120 125 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynadv_ubs.F90
r12377 r13151 34 34 !! * Substitutions 35 35 # include "do_loop_substitute.h90" 36 # include "domzgr_substitute.h90" 36 37 !!---------------------------------------------------------------------- 37 38 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 169 170 DO_2D_00_00 170 171 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 171 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 172 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) & 173 & / e3u(ji,jj,jk,Kmm) 172 174 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 173 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 175 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) & 176 & / e3v(ji,jj,jk,Kmm) 174 177 END_2D 175 178 END DO … … 206 209 END DO 207 210 DO_3D_00_00( 1, jpkm1 ) 208 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 209 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 211 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) & 212 & / e3u(ji,jj,jk,Kmm) 213 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) & 214 & / e3v(ji,jj,jk,Kmm) 210 215 END_3D 211 216 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynatf.F90
r12489 r13151 13 13 !! - ! 2002-10 (C. Talandier, A-M. Treguier) Open boundary cond. 14 14 !! 2.0 ! 2005-11 (V. Garnier) Surface pressure gradient organization 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 15 !! 2.3 ! 2007-07 (D. Storkey) Calls to BDY routines. 16 16 !! 3.2 ! 2009-06 (G. Madec, R.Benshila) re-introduce the vvl option 17 17 !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module … … 22 22 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename dynnxt.F90 -> dynatf.F90. Now just does time filtering. 23 23 !!------------------------------------------------------------------------- 24 24 25 25 !!---------------------------------------------------------------------------------------------- 26 26 !! dyn_atf : apply Asselin time filtering to "now" velocities and vertical scale factors … … 42 42 USE trdken ! trend manager: kinetic energy 43 43 USE isf_oce , ONLY: ln_isf ! ice shelf 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 44 USE isfdynatf , ONLY: isf_dynatf ! ice shelf volume filter correction subroutine 45 45 ! 46 46 USE in_out_manager ! I/O manager … … 59 59 PUBLIC dyn_atf ! routine called by step.F90 60 60 61 #if defined key_qco 62 !!---------------------------------------------------------------------- 63 !! 'key_qco' EMPTY ROUTINE Quasi-Eulerian vertical coordonate 64 !!---------------------------------------------------------------------- 65 CONTAINS 66 67 SUBROUTINE dyn_atf ( kt, Kbb, Kmm, Kaa, puu, pvv, pe3t, pe3u, pe3v ) 68 INTEGER , INTENT(in ) :: kt ! ocean time-step index 69 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! before and after time level indices 70 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered 71 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: pe3t, pe3u, pe3v ! scale factors to be time filtered 72 73 WRITE(*,*) 'dyn_atf: You should not have seen this print! error?', kt 74 END SUBROUTINE dyn_atf 75 76 #else 77 61 78 !! * Substitutions 62 79 # include "do_loop_substitute.h90" 63 80 !!---------------------------------------------------------------------- 64 81 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 65 !! $Id$ 82 !! $Id$ 66 83 !! Software governed by the CeCILL license (see ./LICENSE) 67 84 !!---------------------------------------------------------------------- … … 71 88 !!---------------------------------------------------------------------- 72 89 !! *** ROUTINE dyn_atf *** 73 !! 74 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 90 !! 91 !! ** Purpose : Finalize after horizontal velocity. Apply the boundary 75 92 !! condition on the after velocity and apply the Asselin time 76 93 !! filter to the now fields. … … 79 96 !! estimate (ln_dynspg_ts=T) 80 97 !! 81 !! * Apply lateral boundary conditions on after velocity 98 !! * Apply lateral boundary conditions on after velocity 82 99 !! at the local domain boundaries through lbc_lnk call, 83 100 !! at the one-way open boundaries (ln_bdy=T), … … 86 103 !! * Apply the Asselin time filter to the now fields 87 104 !! arrays to start the next time step: 88 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 105 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 89 106 !! + rn_atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 90 107 !! Note that with flux form advection and non linear free surface, … … 92 109 !! As a result, dyn_atf MUST be called after tra_atf. 93 110 !! 94 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 111 !! ** Action : puu(Kmm),pvv(Kmm) filtered now horizontal velocity 95 112 !!---------------------------------------------------------------------- 96 113 INTEGER , INTENT(in ) :: kt ! ocean time-step index … … 103 120 REAL(wp) :: zve3a, zve3n, zve3b, z1_2dt ! - - 104 121 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve, zwfld 105 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 122 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3t_f, ze3u_f, ze3v_f, zua, zva 106 123 !!---------------------------------------------------------------------- 107 124 ! … … 131 148 ! 132 149 IF( .NOT.ln_bt_fw ) THEN 133 ! Remove advective velocity from "now velocities" 134 ! prior to asselin filtering 135 ! In the forward case, this is done below after asselin filtering 136 ! so that asselin contribution is removed at the same time 150 ! Remove advective velocity from "now velocities" 151 ! prior to asselin filtering 152 ! In the forward case, this is done below after asselin filtering 153 ! so that asselin contribution is removed at the same time 137 154 DO jk = 1, jpkm1 138 155 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 139 156 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 140 END DO 157 END DO 141 158 ENDIF 142 159 ENDIF 143 160 144 161 ! Update after velocity on domain lateral boundaries 145 ! -------------------------------------------------- 162 ! -------------------------------------------------- 146 163 # if defined key_agrif 147 164 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries … … 198 215 zwfld(:,:) = emp_b(:,:) - emp(:,:) 199 216 IF ( ln_rnf ) zwfld(:,:) = zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) ) 217 200 218 DO jk = 1, jpkm1 201 219 ze3t_f(:,:,jk) = ze3t_f(:,:,jk) - zcoef * zwfld(:,:) * tmask(:,:,jk) & 202 & * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 220 & * pe3t(:,:,jk,Kmm) / ( ht(:,:) + 1._wp - ssmask(:,:) ) 203 221 END DO 204 222 ! … … 237 255 pvv(ji,jj,jk,Kmm) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 238 256 END_3D 239 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 257 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) 240 258 pe3v(:,:,1:jpkm1,Kmm) = ze3v_f(:,:,1:jpkm1) 241 259 ! … … 248 266 IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN 249 267 ! Revert filtered "now" velocities to time split estimate 250 ! Doing it here also means that asselin filter contribution is removed 268 ! Doing it here also means that asselin filter contribution is removed 251 269 zue(:,:) = pe3u(:,:,1,Kmm) * puu(:,:,1,Kmm) * umask(:,:,1) 252 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 270 zve(:,:) = pe3v(:,:,1,Kmm) * pvv(:,:,1,Kmm) * vmask(:,:,1) 253 271 DO jk = 2, jpkm1 254 272 zue(:,:) = zue(:,:) + pe3u(:,:,jk,Kmm) * puu(:,:,jk,Kmm) * umask(:,:,jk) 255 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 273 zve(:,:) = zve(:,:) + pe3v(:,:,jk,Kmm) * pvv(:,:,jk,Kmm) * vmask(:,:,jk) 256 274 END DO 257 275 DO jk = 1, jpkm1 … … 305 323 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' nxt - puu(:,:,:,Kaa): ', mask1=umask, & 306 324 & tab3d_2=pvv(:,:,:,Kaa), clinfo2=' pvv(:,:,:,Kaa): ' , mask2=vmask ) 307 ! 325 ! 308 326 IF( ln_dynspg_ts ) DEALLOCATE( zue, zve ) 309 327 IF( l_trddyn ) DEALLOCATE( zua, zva ) … … 312 330 END SUBROUTINE dyn_atf 313 331 332 #endif 333 314 334 !!========================================================================= 315 335 END MODULE dynatf -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynhpg.F90
r12377 r13151 43 43 USE in_out_manager ! I/O manager 44 44 USE prtctl ! Print control 45 USE lbclnk ! lateral boundary condition 45 USE lbclnk ! lateral boundary condition 46 46 USE lib_mpp ! MPP library 47 47 USE eosbn2 ! compute density … … 76 76 !! * Substitutions 77 77 # include "do_loop_substitute.h90" 78 # include "domzgr_substitute.h90" 79 78 80 !!---------------------------------------------------------------------- 79 81 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 204 206 ! 205 207 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 206 ! 208 ! 207 209 IF(lwp) THEN 208 210 WRITE(numout,*) … … 217 219 WRITE(numout,*) 218 220 ENDIF 219 ! 221 ! 220 222 END SUBROUTINE dyn_hpg_init 221 223 … … 427 429 zcpx(ji,jj) = 0._wp 428 430 END IF 429 431 430 432 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 431 433 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 452 454 DO_2D_00_00 453 455 ! hydrostatic pressure gradient along s-surfaces 454 zhpi(ji,jj,1) = zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 455 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e1u(ji,jj) 456 zhpj(ji,jj,1) = zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 457 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) * r1_e2v(ji,jj) 456 zhpi(ji,jj,1) = & 457 & zcoef0 * ( e3w(ji+1,jj ,1,Kmm) * ( znad + rhd(ji+1,jj ,1) ) & 458 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 459 & * r1_e1u(ji,jj) 460 zhpj(ji,jj,1) = & 461 & zcoef0 * ( e3w(ji ,jj+1,1,Kmm) * ( znad + rhd(ji ,jj+1,1) ) & 462 & - e3w(ji ,jj ,1,Kmm) * ( znad + rhd(ji ,jj ,1) ) ) & 463 & * r1_e2v(ji,jj) 458 464 ! s-coordinate pressure gradient correction 459 465 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 464 470 IF( ln_wd_il ) THEN 465 471 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 466 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 472 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 467 473 zuap = zuap * zcpx(ji,jj) 468 474 zvap = zvap * zcpy(ji,jj) … … 478 484 ! hydrostatic pressure gradient along s-surfaces 479 485 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 480 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad )&481 & 486 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 487 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 482 488 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 483 & 484 & 489 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 490 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 485 491 ! s-coordinate pressure gradient correction 486 492 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & … … 491 497 IF( ln_wd_il ) THEN 492 498 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 493 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 499 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 494 500 zuap = zuap * zcpx(ji,jj) 495 501 zvap = zvap * zcpy(ji,jj) … … 522 528 !! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj 523 529 !! iceload is added 524 !! 530 !! 525 531 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend 526 532 !!---------------------------------------------------------------------- … … 540 546 znad=1._wp ! To use density and not density anomaly 541 547 ! 542 ! ! iniitialised to 0. zhpi zhpi 548 ! ! iniitialised to 0. zhpi zhpi 543 549 zhpi(:,:,:) = 0._wp ; zhpj(:,:,:) = 0._wp 544 550 … … 554 560 CALL eos( zts_top, risfdep, zrhdtop_oce ) 555 561 556 !================================================================================== 557 !===== Compute surface value ===================================================== 562 !================================================================================== 563 !===== Compute surface value ===================================================== 558 564 !================================================================================== 559 565 DO_2D_00_00 … … 567 573 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 568 574 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 569 & + ( risfload(ji+1,jj) - risfload(ji,jj)) ) 575 & + ( risfload(ji+1,jj) - risfload(ji,jj)) ) 570 576 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w(ji,jj+1,iktp1j,Kmm) & 571 577 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 572 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 578 & - 0.5_wp * e3w(ji,jj,ikt,Kmm) & 573 579 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 574 & + ( risfload(ji,jj+1) - risfload(ji,jj)) ) 580 & + ( risfload(ji,jj+1) - risfload(ji,jj)) ) 575 581 ! s-coordinate pressure gradient correction (=0 if z coordinate) 576 582 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 582 588 pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 583 589 END_2D 584 !================================================================================== 585 !===== Compute interior value ===================================================== 590 !================================================================================== 591 !===== Compute interior value ===================================================== 586 592 !================================================================================== 587 593 ! interior value (2=<jk=<jpkm1) … … 589 595 ! hydrostatic pressure gradient along s-surfaces 590 596 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 591 & * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 592 & - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 597 & * ( e3w(ji+1,jj,jk,Kmm) & 598 & * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 599 & - e3w(ji ,jj,jk,Kmm) & 600 & * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 593 601 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 594 & * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 595 & - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 602 & * ( e3w(ji,jj+1,jk,Kmm) & 603 & * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 604 & - e3w(ji,jj ,jk,Kmm) & 605 & * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 596 606 ! s-coordinate pressure gradient correction 597 607 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & … … 650 660 zcpx(ji,jj) = 0._wp 651 661 END IF 652 662 653 663 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 654 664 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 771 781 !------------------------------------------------------------- 772 782 773 !!bug gm : e3w-gde3w = 0.5*e3w .... and gde3w(2)-gde3w(1)=e3w(2) .... to be verified774 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be783 !!bug gm : e3w-gde3w(:,:,:) = 0.5*e3w .... and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) .... to be verified 784 ! true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be 775 785 776 786 DO_2D_00_00 … … 825 835 IF( ln_wd_il ) THEN 826 836 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) 827 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 837 zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj) 828 838 ENDIF 829 839 ! add to the general momentum trend … … 845 855 IF( ln_wd_il ) THEN 846 856 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj) 847 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 857 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj) 848 858 ENDIF 849 859 ! add to the general momentum trend … … 916 926 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 917 927 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 918 928 919 929 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 920 930 ELSE 921 931 zcpx(ji,jj) = 0._wp 922 932 END IF 923 933 924 934 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 925 935 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & … … 1002 1012 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1003 1013 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * & 1004 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1014 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1005 1015 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * & 1006 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1016 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1007 1017 !!gm not this: 1008 1018 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1009 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1019 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1010 1020 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1011 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1021 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1012 1022 END_2D 1013 1023 … … 1015 1025 1016 1026 DO_2D_00_00 1017 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1027 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) * znad) 1018 1028 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) * znad) 1019 1029 END_2D … … 1098 1108 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1099 1109 ENDIF 1100 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1110 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) 1101 1111 ENDIF 1102 1112 … … 1154 1164 ENDIF 1155 1165 IF( ln_wd_il ) THEN 1156 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1157 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1166 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1167 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1158 1168 ENDIF 1159 1169 … … 1189 1199 !!---------------------------------------------------------------------- 1190 1200 ! 1191 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1201 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1192 1202 jpi = size(fsp,1) 1193 1203 jpj = size(fsp,2) … … 1359 1369 !!====================================================================== 1360 1370 END MODULE dynhpg 1361 -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynldf_iso.F90
r12377 r13151 22 22 USE ldftra ! lateral physics: eddy diffusivity 23 23 USE zdf_oce ! ocean vertical physics 24 USE ldfslp ! iso-neutral slopes 24 USE ldfslp ! iso-neutral slopes 25 25 ! 26 26 USE in_out_manager ! I/O manager … … 36 36 37 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akzu, akzv !: vertical component of rotated lateral viscosity 38 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso) 38 39 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso) 40 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - - 41 41 42 42 !! * Substitutions 43 43 # include "do_loop_substitute.h90" 44 # include "domzgr_substitute.h90" 44 45 !!---------------------------------------------------------------------- 45 46 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 53 54 !! *** ROUTINE dyn_ldf_iso_alloc *** 54 55 !!---------------------------------------------------------------------- 55 ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , & 56 ALLOCATE( akzu(jpi,jpj,jpk) , zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , & 56 57 & akzv(jpi,jpj,jpk) , zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc ) 57 58 ! … … 63 64 !!---------------------------------------------------------------------- 64 65 !! *** ROUTINE dyn_ldf_iso *** 65 !! 66 !! 66 67 !! ** Purpose : Compute the before trend of the rotated laplacian 67 68 !! operator of lateral momentum diffusion except the diagonal … … 137 138 ! 138 139 ENDIF 139 140 140 141 zaht_0 = 0.5_wp * rn_Ud * rn_Ld ! aht_0 from namtra_ldf = zaht_max 141 142 142 143 ! ! =============== 143 144 DO jk = 1, jpkm1 ! Horizontal slab … … 161 162 162 163 ! -----f----- 163 ! Horizontal fluxes on U | 164 ! Horizontal fluxes on U | 164 165 ! --------------------=== t u t 165 ! | 166 ! | 166 167 ! i-flux at t-point -----f----- 167 168 168 169 IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 169 170 DO_2D_00_01 170 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * MIN( e3u(ji,jj,jk,Kmm), e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 171 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) & 172 & * MIN( e3u(ji ,jj,jk,Kmm), & 173 & e3u(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) 171 174 172 175 zmskt = 1._wp / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) & … … 181 184 ELSE ! other coordinate system (zco or sco) : e3t 182 185 DO_2D_00_01 183 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 186 zabe1 = ( ahmt(ji,jj,jk)+rn_ahm_b ) & 187 & * e2t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e1t(ji,jj) 184 188 185 189 zmskt = 1._wp / MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk+1) & … … 196 200 ! j-flux at f-point 197 201 DO_2D_10_10 198 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 202 zabe2 = ( ahmf(ji,jj,jk) + rn_ahm_b ) & 203 & * e1f(ji,jj) * e3f(ji,jj,jk) * r1_e2f(ji,jj) 199 204 200 205 zmskf = 1._wp / MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) & … … 215 220 216 221 DO_2D_00_10 217 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 222 zabe1 = ( ahmf(ji,jj,jk) + rn_ahm_b ) & 223 & * e2f(ji,jj) * e3f(ji,jj,jk) * r1_e1f(ji,jj) 218 224 219 225 zmskf = 1._wp / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) & … … 230 236 IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u) 231 237 DO_2D_01_10 232 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * MIN( e3v(ji,jj,jk,Kmm), e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 238 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) & 239 & * MIN( e3v(ji,jj ,jk,Kmm), & 240 & e3v(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) 233 241 234 242 zmskt = 1._wp / MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 243 251 ELSE ! other coordinate system (zco or sco) : e3t 244 252 DO_2D_01_10 245 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 253 zabe2 = ( ahmt(ji,jj,jk)+rn_ahm_b ) & 254 & * e1t(ji,jj) * e3t(ji,jj,jk,Kmm) * r1_e2t(ji,jj) 246 255 247 256 zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) & … … 261 270 DO_2D_00_00 262 271 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( ziut(ji+1,jj) - ziut(ji,jj ) & 263 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 272 & + zjuf(ji ,jj) - zjuf(ji,jj-1) ) * r1_e1e2u(ji,jj) & 273 & / e3u(ji,jj,jk,Kmm) 264 274 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zivf(ji,jj ) - zivf(ji-1,jj) & 265 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 275 & + zjvt(ji,jj+1) - zjvt(ji,jj ) ) * r1_e1e2v(ji,jj) & 276 & / e3v(ji,jj,jk,Kmm) 266 277 END_2D 267 278 ! ! =============== … … 278 289 ! ! =============== 279 290 280 291 281 292 ! I. vertical trends associated with the lateral mixing 282 293 ! ===================================================== … … 375 386 DO jk = 1, jpkm1 376 387 DO ji = 2, jpim1 377 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 378 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 388 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + ( zfuw(ji,jk) - zfuw(ji,jk+1) ) * r1_e1e2u(ji,jj) & 389 & / e3u(ji,jj,jk,Kmm) 390 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + ( zfvw(ji,jk) - zfvw(ji,jk+1) ) * r1_e1e2v(ji,jj) & 391 & / e3v(ji,jj,jk,Kmm) 379 392 END DO 380 393 END DO -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynldf_lap_blp.F90
r12377 r13151 14 14 USE dom_oce ! ocean space and time domain 15 15 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 16 USE ldfslp ! iso-neutral slopes 16 USE ldfslp ! iso-neutral slopes 17 17 USE zdf_oce ! ocean vertical physics 18 18 ! … … 28 28 !! * Substitutions 29 29 # include "do_loop_substitute.h90" 30 # include "domzgr_substitute.h90" 30 31 !!---------------------------------------------------------------------- 31 32 !! NEMO/OCE 4.0 , NEMO Consortium (2018) 32 !! $Id$ 33 !! $Id$ 33 34 !! Software governed by the CeCILL license (see ./LICENSE) 34 35 !!---------------------------------------------------------------------- … … 38 39 !!---------------------------------------------------------------------- 39 40 !! *** ROUTINE dyn_ldf_lap *** 40 !! 41 !! ** Purpose : Compute the before horizontal momentum diffusive 41 !! 42 !! ** Purpose : Compute the before horizontal momentum diffusive 42 43 !! trend and add it to the general trend of momentum equation. 43 44 !! 44 !! ** Method : The Laplacian operator apply on horizontal velocity is 45 !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 45 !! ** Method : The Laplacian operator apply on horizontal velocity is 46 !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) 46 47 !! 47 48 !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. … … 76 77 !!gm open question here : e3f at before or now ? probably now... 77 78 !!gm note that ahmf has already been multiplied by fmask 78 zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 79 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 80 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 79 zcur(ji-1,jj-1) = & 80 & ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & 81 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 82 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 81 83 ! ! ahm * div (computed from 2 to jpi/jpj) 82 84 !!gm note that ahmt has already been multiplied by tmask 83 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & 84 & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & 85 & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) 85 zdiv(ji,jj) = & 86 & ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & 87 & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) & 88 & - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & 89 & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) & 90 & - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) 86 91 END_2D 87 92 ! 88 93 DO_2D_00_00 89 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( & 90 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 91 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 94 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * ( & 95 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) & 96 & / e3u(ji,jj,jk,Kmm) & 97 & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) 92 98 ! 93 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * ( & 94 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) & 95 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 99 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * ( & 100 & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) & 101 & / e3v(ji,jj,jk,Kmm) & 102 & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) 96 103 END_2D 97 104 ! ! =============== … … 105 112 !!---------------------------------------------------------------------- 106 113 !! *** ROUTINE dyn_ldf_blp *** 107 !! 108 !! ** Purpose : Compute the before lateral momentum viscous trend 114 !! 115 !! ** Purpose : Compute the before lateral momentum viscous trend 109 116 !! and add it to the general trend of momentum equation. 110 117 !! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynspg_ts.F90
r12489 r13151 87 87 !! * Substitutions 88 88 # include "do_loop_substitute.h90" 89 # include "domzgr_substitute.h90" 89 90 !!---------------------------------------------------------------------- 90 91 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 161 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 162 163 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 164 REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 163 165 ! 164 166 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 227 229 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 228 230 ! ! --------------------------- ! 229 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 230 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 231 DO jk = 1 , jpk 232 ze3u(:,:,jk) = e3u(:,:,jk,Kmm) 233 ze3v(:,:,jk) = e3v(:,:,jk,Kmm) 234 END DO 235 ! 236 zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 237 zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 231 238 ! 232 239 ! … … 250 257 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 251 258 ! 252 CALL dyn_cor_2d( h u(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in259 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 253 260 & zu_trd, zv_trd ) ! ==>> out 254 261 ! … … 567 574 ! at each time step. We however keep them constant here for optimization. 568 575 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 569 CALL dyn_cor_2d( zh up2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )576 CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 570 577 ! 571 578 ! Add tidal astronomical forcing if defined … … 1088 1095 ! 1089 1096 SELECT CASE( nvor_scheme ) 1090 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)1097 CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme 1091 1098 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1092 1099 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) … … 1115 1122 END_2D 1116 1123 ! 1117 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)1124 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme 1118 1125 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1119 1126 DO_2D_01_01 … … 1179 1186 1180 1187 1181 SUBROUTINE dyn_cor_2d( ph u, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )1188 SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd ) 1182 1189 !!--------------------------------------------------------------------- 1183 1190 !! *** ROUTINE dyn_cor_2d *** … … 1187 1194 INTEGER :: ji ,jj ! dummy loop indices 1188 1195 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1189 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ph u, phv, punb, pvnb, zhU, zhV1196 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV 1190 1197 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1191 1198 !!---------------------------------------------------------------------- … … 1196 1203 z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1197 1204 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1198 & * ( e1e2t(ji+1,jj)* ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &1199 & + e1e2t(ji ,jj)* ht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )1205 & * ( e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) & 1206 & + e1e2t(ji ,jj)*pht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) ) 1200 1207 ! 1201 1208 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1202 & * ( e1e2t(ji,jj+1)* ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &1203 & + e1e2t(ji,jj )* ht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )1209 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1210 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1204 1211 END_2D 1205 1212 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynvor.F90
r12377 r13151 15 15 !! 3.2 ! 2009-04 (R. Benshila) vvl: correction of een scheme 16 16 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 17 !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 17 !! 3.7 ! 2014-04 (G. Madec) trend simplification: suppress jpdyn_trd_dat vorticity 18 18 !! - ! 2014-06 (G. Madec) suppression of velocity curl from in-core memory 19 19 !! - ! 2016-12 (G. Madec, E. Clementi) add Stokes-Coriolis trends (ln_stcor=T) … … 70 70 INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme 71 71 72 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 72 INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity 73 73 ! ! associated indices: 74 74 INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary) … … 79 79 80 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 - - - - 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1v_2 ! = dj(e1v)/2 - - - - 82 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) - - - - 84 83 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: dj_e1u_2e1e2f ! = dj(e1v)/(2*e1e2f) - - - - 84 85 85 REAL(wp) :: r1_4 = 0.250_wp ! =1/4 86 86 REAL(wp) :: r1_8 = 0.125_wp ! =1/8 87 87 REAL(wp) :: r1_12 = 1._wp / 12._wp ! 1/12 88 88 89 89 !! * Substitutions 90 90 # include "do_loop_substitute.h90" 91 # include "domzgr_substitute.h90" 92 91 93 !!---------------------------------------------------------------------- 92 94 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 103 105 !! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend 104 106 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative 105 !! and planetary vorticity trends) and send them to trd_dyn 107 !! and planetary vorticity trends) and send them to trd_dyn 106 108 !! for futher diagnostics (l_trddyn=T) 107 109 !!---------------------------------------------------------------------- … … 193 195 !! *** ROUTINE vor_enT *** 194 196 !! 195 !! ** Purpose : Compute the now total vorticity trend and add it to 197 !! ** Purpose : Compute the now total vorticity trend and add it to 196 198 !! the general trend of the momentum equation. 197 199 !! 198 !! ** Method : Trend evaluated using now fields (centered in time) 200 !! ** Method : Trend evaluated using now fields (centered in time) 199 201 !! and t-point evaluation of vorticity (planetary and relative). 200 202 !! conserves the horizontal kinetic energy. 201 !! The general trend of momentum is increased due to the vorticity 203 !! The general trend of momentum is increased due to the vorticity 202 204 !! term which is given by: 203 205 !! voru = 1/bu mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t mj[vn] ] … … 233 235 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 234 236 END_2D 235 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 237 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 236 238 DO_2D_10_10 237 239 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) … … 248 250 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 249 251 END_2D 250 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 252 IF( ln_dynvor_msk ) THEN ! mask/unmask relative vorticity 251 253 DO_2D_10_10 252 254 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) … … 269 271 DO_2D_01_01 270 272 zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 271 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 273 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) & 274 & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 272 275 END_2D 273 276 CASE ( np_MET ) !* metric term 274 277 DO_2D_01_01 275 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 276 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm) 278 zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 279 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) & 280 & * e3t(ji,jj,jk,Kmm) 277 281 END_2D 278 282 CASE ( np_CRV ) !* Coriolis + relative vorticity 279 283 DO_2D_01_01 280 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 281 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 284 zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ,jk) + zwz(ji,jj ,jk) & 285 & + zwz(ji-1,jj-1,jk) + zwz(ji,jj-1,jk) ) ) & 286 & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 282 287 END_2D 283 288 CASE ( np_CME ) !* Coriolis + metric 284 289 DO_2D_01_01 285 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 286 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 287 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) * e3t(ji,jj,jk,Kmm) 290 zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) & 291 & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) & 292 & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) & 293 & * e3t(ji,jj,jk,Kmm) 288 294 END_2D 289 295 CASE DEFAULT ! error … … 298 304 ! 299 305 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) & 300 & * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) & 301 & + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) ) 306 & * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) & 307 & + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) ) 302 308 END_2D 303 309 ! ! =============== … … 311 317 !! *** ROUTINE vor_ene *** 312 318 !! 313 !! ** Purpose : Compute the now total vorticity trend and add it to 319 !! ** Purpose : Compute the now total vorticity trend and add it to 314 320 !! the general trend of the momentum equation. 315 321 !! 316 !! ** Method : Trend evaluated using now fields (centered in time) 322 !! ** Method : Trend evaluated using now fields (centered in time) 317 323 !! and the Sadourny (1975) flux form formulation : conserves the 318 324 !! horizontal kinetic energy. 319 !! The general trend of momentum is increased due to the vorticity 325 !! The general trend of momentum is increased due to the vorticity 320 326 !! term which is given by: 321 327 !! voru = 1/e1u mj-1[ (rvor+f)/e3f mi(e1v*e3v pvv(:,:,:,Kmm)) ] … … 350 356 SELECT CASE( kvor ) !== vorticity considered ==! 351 357 CASE ( np_COR ) !* Coriolis (planetary vorticity) 352 zwz(:,:) = ff_f(:,:) 358 zwz(:,:) = ff_f(:,:) 353 359 CASE ( np_RVO ) !* relative vorticity 354 360 DO_2D_10_10 … … 396 402 zx2 = zwx(ji ,jj) + zwx(ji ,jj+1) 397 403 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1u(ji,jj) * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 398 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 404 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e2v(ji,jj) * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 399 405 END_2D 400 406 ! ! =============== … … 446 452 SELECT CASE( kvor ) !== vorticity considered ==! 447 453 CASE ( np_COR ) !* Coriolis (planetary vorticity) 448 zwz(:,:) = ff_f(:,:) 454 zwz(:,:) = ff_f(:,:) 449 455 CASE ( np_RVO ) !* relative vorticity 450 456 DO_2D_10_10 … … 504 510 !! *** ROUTINE vor_een *** 505 511 !! 506 !! ** Purpose : Compute the now total vorticity trend and add it to 512 !! ** Purpose : Compute the now total vorticity trend and add it to 507 513 !! the general trend of the momentum equation. 508 514 !! 509 !! ** Method : Trend evaluated using now fields (centered in time) 510 !! and the Arakawa and Lamb (1980) flux form formulation : conserves 515 !! ** Method : Trend evaluated using now fields (centered in time) 516 !! and the Arakawa and Lamb (1980) flux form formulation : conserves 511 517 !! both the horizontal kinetic energy and the potential enstrophy 512 518 !! when horizontal divergence is zero (see the NEMO documentation) … … 545 551 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 546 552 DO_2D_10_10 547 ze3f = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 548 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 553 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 554 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 555 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 556 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 549 557 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f 550 558 ELSE ; z1_e3f(ji,jj) = 0._wp … … 553 561 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 554 562 DO_2D_10_10 555 ze3f = ( e3t(ji,jj+1,jk,Kmm)*tmask(ji,jj+1,jk) + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 556 & + e3t(ji,jj ,jk,Kmm)*tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 563 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 564 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 565 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 566 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 557 567 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 558 568 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) … … 644 654 !! *** ROUTINE vor_eeT *** 645 655 !! 646 !! ** Purpose : Compute the now total vorticity trend and add it to 656 !! ** Purpose : Compute the now total vorticity trend and add it to 647 657 !! the general trend of the momentum equation. 648 658 !! 649 !! ** Method : Trend evaluated using now fields (centered in time) 650 !! and the Arakawa and Lamb (1980) vector form formulation using 659 !! ** Method : Trend evaluated using now fields (centered in time) 660 !! and the Arakawa and Lamb (1980) vector form formulation using 651 661 !! a modified version of Arakawa and Lamb (1980) scheme (see vor_een). 652 !! The change consists in 662 !! The change consists in 653 663 !! Add this trend to the general momentum trend (pu_rhs,pv_rhs). 654 664 !! … … 667 677 REAL(wp) :: zua, zva ! local scalars 668 678 REAL(wp) :: zmsk, z1_e3t ! local scalars 669 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy 679 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy 670 680 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 671 681 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz … … 827 837 ! 828 838 IF( ioptio /= 1 ) CALL ctl_stop( ' use ONE and ONLY one vorticity scheme' ) 829 ! 839 ! 830 840 IF(lwp) WRITE(numout,*) ! type of calculated vorticity (set ncor, nrvm, ntot) 831 841 ncor = np_COR ! planetary vorticity … … 836 846 ntot = np_COR ! - - 837 847 CASE( np_VEC_c2 ) 838 IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity' 848 IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity' 839 849 nrvm = np_RVO ! relative vorticity 840 ntot = np_CRV ! relative + planetary vorticity 850 ntot = np_CRV ! relative + planetary vorticity 841 851 CASE( np_FLX_c2 , np_FLX_ubs ) 842 852 IF(lwp) WRITE(numout,*) ' ==>>> flux form dynamics : total vorticity = Coriolis + metric term' … … 863 873 ! 864 874 END SELECT 865 875 866 876 IF(lwp) THEN ! Print the choice 867 877 WRITE(numout,*) … … 873 883 CASE( np_EEN ) ; WRITE(numout,*) ' ==>>> energy and enstrophy conserving scheme (EEN)' 874 884 CASE( np_MIX ) ; WRITE(numout,*) ' ==>>> mixed enstrophy/energy conserving scheme (MIX)' 875 END SELECT 885 END SELECT 876 886 ENDIF 877 887 ! -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynzad.F90
r12377 r13151 29 29 !! * Substitutions 30 30 # include "do_loop_substitute.h90" 31 # include "domzgr_substitute.h90" 31 32 !!---------------------------------------------------------------------- 32 33 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 95 96 ! 96 97 DO_3D_00_00( 1, jpkm1 ) 97 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) 98 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) 98 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) & 99 & / e3u(ji,jj,jk,Kmm) 100 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) & 101 & / e3v(ji,jj,jk,Kmm) 99 102 END_3D 100 103 -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/dynzdf.F90
r12489 r13151 38 38 !! * Substitutions 39 39 # include "do_loop_substitute.h90" 40 # include "domzgr_substitute.h90" 40 41 !!---------------------------------------------------------------------- 41 42 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 55 56 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 57 !! u(after) = u(before) + 2*dt * u(rhs) vector form or linear free surf. 57 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u (after)otherwise58 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after otherwise 58 59 !! - update the after velocity with the implicit vertical mixing. 59 60 !! This requires to solver the following system: 60 !! u(after) = u(after) + 1/e3u (after) dk+1[ mi(avm) / e3uw(after)dk[ua] ]61 !! u(after) = u(after) + 1/e3u_after dk+1[ mi(avm) / e3uw_after dk[ua] ] 61 62 !! with the following surface/top/bottom boundary condition: 62 63 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) … … 112 113 ELSE ! applied on thickness weighted velocity 113 114 DO jk = 1, jpkm1 114 puu(:,:,jk,Kaa) = ( e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) & 115 & + rDt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) ) / e3u(:,:,jk,Kaa) * umask(:,:,jk) 116 pvv(:,:,jk,Kaa) = ( e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) & 117 & + rDt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk) 115 puu(:,:,jk,Kaa) = ( e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) & 116 & + rDt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) ) & 117 & / e3u(:,:,jk,Kaa) * umask(:,:,jk) 118 pvv(:,:,jk,Kaa) = ( e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) & 119 & + rDt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) ) & 120 & / e3v(:,:,jk,Kaa) * vmask(:,:,jk) 118 121 END DO 119 122 ENDIF … … 131 134 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 132 135 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 133 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 134 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 136 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 137 & + r_vvl * e3u(ji,jj,iku,Kaa) 138 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 139 & + r_vvl * e3v(ji,jj,ikv,Kaa) 135 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) / ze3ua 136 141 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) / ze3va … … 140 145 iku = miku(ji,jj) ! top ocean level at u- and v-points 141 146 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 142 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 143 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 147 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 148 & + r_vvl * e3u(ji,jj,iku,Kaa) 149 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 150 & + r_vvl * e3v(ji,jj,ikv,Kaa) 144 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) / ze3ua 145 152 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) / ze3va … … 156 163 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 157 164 DO_3D_00_00( 1, jpkm1 ) 158 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 165 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 166 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 159 167 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 160 168 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) … … 169 177 CASE DEFAULT ! iso-level lateral mixing 170 178 DO_3D_00_00( 1, jpkm1 ) 171 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 172 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 173 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 179 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & ! after scale factor at U-point 180 & + r_vvl * e3u(ji,jj,jk,Kaa) 181 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 182 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 183 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 184 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 174 185 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 175 186 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua … … 181 192 DO_2D_00_00 182 193 zwi(ji,jj,1) = 0._wp 183 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 184 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 194 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 195 & + r_vvl * e3u(ji,jj,1,Kaa) 196 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) & 197 & / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 185 198 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 186 199 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) … … 191 204 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 192 205 DO_3D_00_00( 1, jpkm1 ) 193 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 206 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 207 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 194 208 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 195 209 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) … … 202 216 CASE DEFAULT ! iso-level lateral mixing 203 217 DO_3D_00_00( 1, jpkm1 ) 204 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 205 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 206 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 218 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 219 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 220 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 221 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 222 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 223 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 207 224 zwi(ji,jj,jk) = zzwi 208 225 zws(ji,jj,jk) = zzws … … 226 243 DO_2D_00_00 227 244 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 228 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 245 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 246 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 229 247 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 230 248 END_2D … … 233 251 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 234 252 iku = miku(ji,jj) ! ocean top level at u- and v-points 235 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 253 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 254 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 236 255 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 237 256 END_2D … … 259 278 ! 260 279 DO_2D_00_00 261 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 280 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 281 & + r_vvl * e3u(ji,jj,1,Kaa) 262 282 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 263 283 & / ( ze3ua * rho0 ) * umask(ji,jj,1) … … 282 302 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv) 283 303 DO_3D_00_00( 1, jpkm1 ) 284 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 304 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 305 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 285 306 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 286 307 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) … … 295 316 CASE DEFAULT ! iso-level lateral mixing 296 317 DO_3D_00_00( 1, jpkm1 ) 297 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 298 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 299 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 318 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 319 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 321 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 322 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 323 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 300 324 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 301 325 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va … … 307 331 DO_2D_00_00 308 332 zwi(ji,jj,1) = 0._wp 309 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 310 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 333 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 334 & + r_vvl * e3v(ji,jj,1,Kaa) 335 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) & 336 & / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 311 337 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 312 338 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) … … 317 343 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 318 344 DO_3D_00_00( 1, jpkm1 ) 319 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 345 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 346 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 347 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 321 348 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) … … 328 355 CASE DEFAULT ! iso-level lateral mixing 329 356 DO_3D_00_00( 1, jpkm1 ) 330 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 331 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 332 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 357 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 358 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 359 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 360 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 361 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 362 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 333 363 zwi(ji,jj,jk) = zzwi 334 364 zws(ji,jj,jk) = zzws … … 351 381 DO_2D_00_00 352 382 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 353 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 383 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 384 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 354 385 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 355 386 END_2D … … 357 388 DO_2D_00_00 358 389 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 359 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 390 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 391 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 360 392 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 361 393 END_2D … … 383 415 ! 384 416 DO_2D_00_00 385 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 417 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 418 & + r_vvl * e3v(ji,jj,1,Kaa) 386 419 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 387 420 & / ( ze3va * rho0 ) * vmask(ji,jj,1) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/sshwzv.F90
r12489 r13151 1 MODULE sshwzv 1 MODULE sshwzv 2 2 !!============================================================================== 3 3 !! *** MODULE sshwzv *** … … 5 5 !!============================================================================== 6 6 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 7 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 7 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 8 8 !! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module … … 20 20 USE oce ! ocean dynamics and tracers variables 21 21 USE isf_oce ! ice shelf 22 USE dom_oce ! ocean space and time domain variables 22 USE dom_oce ! ocean space and time domain variables 23 23 USE sbc_oce ! surface boundary condition: ocean 24 24 USE domvvl ! Variable volume … … 31 31 #endif 32 32 ! 33 USE iom 33 USE iom 34 34 USE in_out_manager ! I/O manager 35 35 USE restart ! only for lrst_oce … … 50 50 !! * Substitutions 51 51 # include "do_loop_substitute.h90" 52 # include "domzgr_substitute.h90" 53 52 54 !!---------------------------------------------------------------------- 53 55 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 60 62 !!---------------------------------------------------------------------- 61 63 !! *** ROUTINE ssh_nxt *** 62 !! 64 !! 63 65 !! ** Purpose : compute the after ssh (ssh(Kaa)) 64 66 !! … … 74 76 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index 75 77 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height 76 ! 78 ! 77 79 INTEGER :: jk ! dummy loop index 78 80 REAL(wp) :: zcoef ! local scalar … … 106 108 ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to 107 109 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 108 ! 110 ! 109 111 pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 110 112 ! 111 113 #if defined key_agrif 112 Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt ) 114 Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa 115 CALL agrif_ssh( kt ) 113 116 #endif 114 117 ! … … 129 132 END SUBROUTINE ssh_nxt 130 133 131 132 SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa)134 135 SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww ) 133 136 !!---------------------------------------------------------------------- 134 137 !! *** ROUTINE wzv *** 135 !! 138 !! 136 139 !! ** Purpose : compute the now vertical velocity 137 140 !! 138 !! ** Method : - Using the incompressibility hypothesis, the vertical 139 !! velocity is computed by integrating the horizontal divergence 141 !! ** Method : - Using the incompressibility hypothesis, the vertical 142 !! velocity is computed by integrating the horizontal divergence 140 143 !! from the bottom to the surface minus the scale factor evolution. 141 144 !! The boundary conditions are w=0 at the bottom (no flux) and. … … 147 150 INTEGER , INTENT(in) :: kt ! time step 148 151 INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices 149 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! now vertical velocity152 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm 150 153 ! 151 154 INTEGER :: ji, jj, jk ! dummy loop indices … … 160 163 IF(lwp) WRITE(numout,*) '~~~~~ ' 161 164 ! 162 pww(:,:,jpk) = 0._wp 165 pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 163 166 ENDIF 164 167 ! !------------------------------! … … 166 169 ! !------------------------------! 167 170 ! 168 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 169 ALLOCATE( zhdiv(jpi,jpj,jpk) ) 171 ! !===============================! 172 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==! 173 ! !===============================! 174 ALLOCATE( zhdiv(jpi,jpj,jpk) ) 170 175 ! 171 176 DO jk = 1, jpkm1 … … 181 186 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 182 187 ! computation of w 183 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk) & 184 & + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 188 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 189 & + zhdiv(:,:,jk) & 190 & + r1_Dt * ( e3t(:,:,jk,Kaa) & 191 & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 185 192 END DO 186 193 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 187 DEALLOCATE( zhdiv ) 188 ELSE ! z_star and linear free surface cases 189 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 190 ! computation of w 191 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 192 & + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 194 DEALLOCATE( zhdiv ) 195 ! !=================================! 196 ELSEIF( ln_linssh ) THEN !== linear free surface cases ==! 197 ! !=================================! 198 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 199 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk) 200 END DO 201 ! !==========================================! 202 ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco') 203 ! !==========================================! 204 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 205 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 206 & + r1_Dt * ( e3t(:,:,jk,Kaa) & 207 & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 193 208 END DO 194 209 ENDIF … … 200 215 ENDIF 201 216 ! 202 #if defined key_agrif 203 IF( .NOT. AGRIF_Root() ) THEN 204 IF ((nbondi == 1).OR.(nbondi == 2)) pww(nlci-1 , : ,:) = 0.e0 ! east 205 IF ((nbondi == -1).OR.(nbondi == 2)) pww(2 , : ,:) = 0.e0 ! west 206 IF ((nbondj == 1).OR.(nbondj == 2)) pww(: ,nlcj-1 ,:) = 0.e0 ! north 207 IF ((nbondj == -1).OR.(nbondj == 2)) pww(: ,2 ,:) = 0.e0 ! south 208 ENDIF 209 #endif 217 #if defined key_agrif 218 IF( .NOT. AGRIF_Root() ) THEN 219 IF ((nbondi == 1).OR.(nbondi == 2)) pww(nlci-1 , : ,:) = 0.e0 ! east 220 IF ((nbondi == -1).OR.(nbondi == 2)) pww(2 , : ,:) = 0.e0 ! west 221 IF ((nbondj == 1).OR.(nbondj == 2)) pww(: ,nlcj-1 ,:) = 0.e0 ! north 222 IF ((nbondj == -1).OR.(nbondj == 2)) pww(: ,2 ,:) = 0.e0 ! south 223 ENDIF 224 #endif 210 225 ! 211 226 IF( ln_timing ) CALL timing_stop('wzv') … … 214 229 215 230 216 SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )231 SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh, pssh_f ) 217 232 !!---------------------------------------------------------------------- 218 233 !! *** ROUTINE ssh_atf *** … … 229 244 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 230 245 !!---------------------------------------------------------------------- 231 INTEGER , INTENT(in ) :: kt ! ocean time-step index 232 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices 233 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field 246 INTEGER , INTENT(in ) :: kt ! ocean time-step index 247 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices 248 REAL(wp), DIMENSION(jpi,jpj,jpt) , TARGET, INTENT(inout) :: pssh ! SSH field 249 REAL(wp), DIMENSION(jpi,jpj ), OPTIONAL, TARGET, INTENT( out) :: pssh_f ! filtered SSH field 234 250 ! 235 251 REAL(wp) :: zcoef ! local scalar 252 REAL(wp), POINTER, DIMENSION(:,:) :: zssh ! pointer for filtered SSH 236 253 !!---------------------------------------------------------------------- 237 254 ! … … 245 262 ! !== Euler time-stepping: no filter, just swap ==! 246 263 IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps 264 IF( PRESENT( pssh_f ) ) THEN ; zssh => pssh_f 265 ELSE ; zssh => pssh(:,:,Kmm) 266 ENDIF 247 267 ! ! filtered "now" field 248 268 pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) … … 266 286 END SUBROUTINE ssh_atf 267 287 288 268 289 SUBROUTINE wAimp( kt, Kmm ) 269 290 !!---------------------------------------------------------------------- 270 291 !! *** ROUTINE wAimp *** 271 !! 292 !! 272 293 !! ** Purpose : compute the Courant number and partition vertical velocity 273 294 !! if a proportion needs to be treated implicitly 274 295 !! 275 !! ** Method : - 296 !! ** Method : - 276 297 !! 277 298 !! ** action : ww : now vertical velocity (to be handled explicitly) … … 279 300 !! 280 301 !! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent 281 !! implicit scheme for vertical advection in oceanic modeling. 302 !! implicit scheme for vertical advection in oceanic modeling. 282 303 !! Ocean Modelling, 91, 38-69. 283 304 !!---------------------------------------------------------------------- … … 306 327 DO_3D_00_00( 1, jpkm1 ) 307 328 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 308 ! 2*rn_Dt and not rDt (for restartability) 309 Cu_adv(ji,jj,jk) = 2._wp * rn_Dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 310 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 311 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) & 312 & * r1_e1e2t(ji,jj) & 313 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - & 314 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) & 315 & * r1_e1e2t(ji,jj) & 316 & ) * z1_e3t 329 ! 2*rdt and not r2dt (for restartability) 330 Cu_adv(ji,jj,jk) = 2._wp * rDt * & 331 & ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 332 & + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) & 333 & * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 334 & MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) & 335 & * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) & 336 & * r1_e1e2t(ji ,jj) & 337 & + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) & 338 & * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - & 339 & MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) & 340 & * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) & 341 & * r1_e1e2t(ji,jj ) & 342 & ) * z1_e3t 317 343 END_3D 318 344 ELSE 319 345 DO_3D_00_00( 1, jpkm1 ) 320 346 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 321 ! 2*rn_Dt and not rDt (for restartability) 322 Cu_adv(ji,jj,jk) = 2._wp * rn_Dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 323 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 324 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 325 & * r1_e1e2t(ji,jj) & 326 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 327 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 328 & * r1_e1e2t(ji,jj) & 329 & ) * z1_e3t 347 ! 2*rdt and not r2dt (for restartability) 348 Cu_adv(ji,jj,jk) = 2._wp * rDt * & 349 & ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 350 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 351 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 352 & * r1_e1e2t(ji,jj) & 353 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 354 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 355 & * r1_e1e2t(ji,jj) & 356 & ) * z1_e3t 330 357 END_3D 331 358 ENDIF … … 339 366 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) ) 340 367 ! alt: 341 ! IF ( ww(ji,jj,jk) > 0._wp ) THEN 342 ! zCu = Cu_adv(ji,jj,jk) 368 ! IF ( ww(ji,jj,jk) > 0._wp ) THEN 369 ! zCu = Cu_adv(ji,jj,jk) 343 370 ! ELSE 344 371 ! zCu = Cu_adv(ji,jj,jk-1) 345 ! ENDIF 372 ! ENDIF 346 373 ! 347 374 IF( zCu <= Cu_min ) THEN !<-- Fully explicit … … 360 387 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl 361 388 END_3D 362 Cu_adv(:,:,1) = 0._wp 389 Cu_adv(:,:,1) = 0._wp 363 390 ELSE 364 391 ! Fully explicit everywhere … … 366 393 wi (:,:,:) = 0._wp 367 394 ENDIF 368 CALL iom_put("wimp",wi) 395 CALL iom_put("wimp",wi) 369 396 CALL iom_put("wi_cff",Cu_adv) 370 397 CALL iom_put("wexp",ww) -
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/OCE/DYN/wet_dry.F90
r12489 r13151 33 33 !! * Substitutions 34 34 # include "do_loop_substitute.h90" 35 # include "domzgr_substitute.h90" 35 36 !!---------------------------------------------------------------------- 36 37 !! critical depths,filters, limiters,and masks for Wetting and Drying
Note: See TracChangeset
for help on using the changeset viewer.
