Changeset 4990 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Timestamp:
- 2014-12-15T17:42:49+01:00 (9 years ago)
- Location:
- trunk/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 21 edited
-
divcur.F90 (modified) (4 diffs)
-
dynadv.F90 (modified) (6 diffs)
-
dynadv_cen2.F90 (modified) (3 diffs)
-
dynadv_ubs.F90 (modified) (3 diffs)
-
dynbfr.F90 (modified) (5 diffs)
-
dynhpg.F90 (modified) (18 diffs)
-
dynkeg.F90 (modified) (3 diffs)
-
dynldf.F90 (modified) (5 diffs)
-
dynldf_bilap.F90 (modified) (2 diffs)
-
dynldf_bilapg.F90 (modified) (3 diffs)
-
dynldf_iso.F90 (modified) (1 diff)
-
dynldf_lap.F90 (modified) (2 diffs)
-
dynnxt.F90 (modified) (13 diffs)
-
dynspg.F90 (modified) (4 diffs)
-
dynspg_exp.F90 (modified) (1 diff)
-
dynspg_flt.F90 (modified) (9 diffs)
-
dynvor.F90 (modified) (16 diffs)
-
dynzad.F90 (modified) (6 diffs)
-
dynzdf.F90 (modified) (2 diffs)
-
dynzdf_imp.F90 (modified) (12 diffs)
-
sshwzv.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/divcur.F90
r4147 r4990 25 25 USE oce ! ocean dynamics and tracers 26 26 USE dom_oce ! ocean space and time domain 27 USE sbc_oce, ONLY : ln_rnf 27 USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean 28 28 USE sbcrnf ! river runoff 29 USE sbcisf ! ice shelf 29 30 USE cla ! cross land advection (cla_div routine) 30 31 USE in_out_manager ! I/O manager … … 66 67 !! - compute the now divergence given by : 67 68 !! hdivn = 1/(e1t*e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] ) 68 !! correct hdiv with runoff inflow (div_rnf) and cross land flow (div_cla) 69 !! correct hdiv with runoff inflow (div_rnf), ice shelf melting (div_isf) 70 !! and cross land flow (div_cla) 69 71 !! II. vorticity : 70 72 !! - save the curl computed at the previous time-step … … 225 227 ! ! =============== 226 228 227 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 228 IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 ) CALL cla_div ( kt ) ! Cross Land Advection (Update Hor. divergence) 229 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 230 IF( ln_divisf .AND. (nn_isf /= 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field) 231 IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (Update Hor. divergence) 229 232 230 233 ! 4. Lateral boundary conditions on hdivn and rotn … … 323 326 ! ! =============== 324 327 325 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 328 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) ! runoffs (update hdivn field) 329 IF( ln_divisf .AND. (nn_isf .GT. 0) ) CALL sbc_isf_div( hdivn ) ! ice shelf (update hdivn field) 326 330 IF( nn_cla == 1 ) CALL cla_div ( kt ) ! Cross Land Advection (update hdivn field) 327 331 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv.F90
r4624 r4990 30 30 LOGICAL, PUBLIC :: ln_dynadv_cen2 !: flux form - 2nd order centered scheme flag 31 31 LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag 32 LOGICAL, PUBLIC :: ln_dynzad_zts !: vertical advection with sub-timestepping (requires vector form) 32 33 33 34 INTEGER :: nadv ! choice of the formulation and scheme for the advection … … 64 65 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy 65 66 CALL dyn_zad ( kt ) ! vector form : vertical advection 66 CASE ( 1 ) 67 CASE ( 1 ) 68 CALL dyn_keg ( kt ) ! vector form : horizontal gradient of kinetic energy 69 CALL dyn_zad_zts ( kt ) ! vector form : vertical advection with sub-timestepping 70 CASE ( 2 ) 67 71 CALL dyn_adv_cen2( kt ) ! 2nd order centered scheme 68 CASE ( 2)72 CASE ( 3 ) 69 73 CALL dyn_adv_ubs ( kt ) ! 3rd order UBS scheme 70 74 ! … … 91 95 INTEGER :: ios ! Local integer output status for namelist read 92 96 !! 93 NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs 97 NAMELIST/namdyn_adv/ ln_dynadv_vec, ln_dynadv_cen2 , ln_dynadv_ubs, ln_dynzad_zts 94 98 !!---------------------------------------------------------------------- 95 99 … … 111 115 WRITE(numout,*) ' 2nd order centred advection scheme ln_dynadv_cen2 = ', ln_dynadv_cen2 112 116 WRITE(numout,*) ' 3rd order UBS advection scheme ln_dynadv_ubs = ', ln_dynadv_ubs 117 WRITE(numout,*) ' Sub timestepping of vertical advection ln_dynzad_zts = ', ln_dynzad_zts 113 118 ENDIF 114 119 … … 120 125 121 126 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE advection scheme in namelist namdyn_adv' ) 127 IF( ln_dynzad_zts .AND. .NOT. ln_dynadv_vec ) & 128 CALL ctl_stop( 'Sub timestepping of vertical advection requires vector form; set ln_dynadv_vec = .TRUE.' ) 122 129 123 130 ! ! Set nadv 124 131 IF( ln_dynadv_vec ) nadv = 0 125 IF( ln_dynadv_cen2 ) nadv = 1 126 IF( ln_dynadv_ubs ) nadv = 2 132 IF( ln_dynzad_zts ) nadv = 1 133 IF( ln_dynadv_cen2 ) nadv = 2 134 IF( ln_dynadv_ubs ) nadv = 3 127 135 IF( lk_esopa ) nadv = -1 128 136 … … 130 138 WRITE(numout,*) 131 139 IF( nadv == 0 ) WRITE(numout,*) ' vector form : keg + zad + vor is used' 132 IF( nadv == 1 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 133 IF( nadv == 2 ) WRITE(numout,*) ' flux form : UBS scheme is used' 140 IF( nadv == 1 ) WRITE(numout,*) ' vector form : keg + zad_zts + vor is used' 141 IF( nadv == 2 ) WRITE(numout,*) ' flux form : 2nd order scheme is used' 142 IF( nadv == 3 ) WRITE(numout,*) ' flux form : UBS scheme is used' 134 143 IF( nadv == -1 ) WRITE(numout,*) ' esopa test: use all advection formulation' 135 144 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_cen2.F90
r3294 r4990 15 15 USE oce ! ocean dynamics and tracers 16 16 USE dom_oce ! ocean space and time domain 17 USE trdmod_oce ! ocean variables trends 18 USE trdmod ! ocean dynamics trends 17 USE trd_oce ! trends: ocean variables 18 USE trddyn ! trend manager: dynamics 19 ! 19 20 USE in_out_manager ! I/O manager 20 21 USE lib_mpp ! MPP library 21 22 USE prtctl ! Print control 22 USE wrk_nemo 23 USE timing 23 USE wrk_nemo ! Memory Allocation 24 USE timing ! Timing 24 25 25 26 IMPLICIT NONE … … 103 104 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 104 105 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 105 CALL trd_ mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )106 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 106 107 zfu_t(:,:,:) = ua(:,:,:) 107 108 zfv_t(:,:,:) = va(:,:,:) … … 153 154 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 154 155 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 155 CALL trd_ mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )156 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 156 157 ENDIF 157 158 ! ! Control print -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r4153 r4990 16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain 18 USE trdmod ! ocean dynamics trends 19 USE trdmod_oce ! ocean variables trends 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 ! 20 21 USE in_out_manager ! I/O manager 21 22 USE prtctl ! Print control 22 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 24 USE lib_mpp ! MPP library 24 USE wrk_nemo 25 USE timing 25 USE wrk_nemo ! Memory Allocation 26 USE timing ! Timing 26 27 27 28 IMPLICIT NONE … … 196 197 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 197 198 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 198 CALL trd_ mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt )199 CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt ) 199 200 zfu_t(:,:,:) = ua(:,:,:) 200 201 zfv_t(:,:,:) = va(:,:,:) … … 245 246 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 246 247 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 247 CALL trd_ mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt )248 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 248 249 ENDIF 249 250 ! ! Control print -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r3294 r4990 10 10 11 11 !!---------------------------------------------------------------------- 12 !! dyn_bfr : Update the momentum trend with the bottom friction contribution12 !! dyn_bfr : Update the momentum trend with the bottom friction contribution 13 13 !!---------------------------------------------------------------------- 14 USE oce 15 USE dom_oce 16 USE zdf_oce 17 USE zdfbfr 18 USE trd mod ! ocean active dynamics and tracers trends19 USE trd mod_oce ! ocean variables trends20 USE in_out_manager 21 USE prtctl 22 USE timing 23 USE wrk_nemo 14 USE oce ! ocean dynamics and tracers variables 15 USE dom_oce ! ocean space and time domain variables 16 USE zdf_oce ! ocean vertical physics variables 17 USE zdfbfr ! ocean bottom friction variables 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 USE in_out_manager ! I/O manager 21 USE prtctl ! Print control 22 USE timing ! Timing 23 USE wrk_nemo ! Memory Allocation 24 24 25 25 IMPLICIT NONE 26 26 PRIVATE 27 27 28 PUBLIC dyn_bfr 28 PUBLIC dyn_bfr ! routine called by step.F90 29 29 30 30 !! * Substitutions … … 57 57 IF( nn_timing == 1 ) CALL timing_start('dyn_bfr') 58 58 ! 59 !!gm issue: better to put the logical in step to control the call of zdf_bfr 60 !! ==> change the logical from ln_bfrimp to ln_bfr_exp !! 59 61 IF( .NOT.ln_bfrimp) THEN ! only for explicit bottom friction form 60 62 ! implicit bfr is implemented in dynzdf_imp 61 63 64 !!gm bug : time step is only rdt (not 2 rdt if euler start !) 62 65 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 63 66 … … 69 72 70 73 71 # if defined key_vectopt_loop72 DO jj = 1, 173 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)74 # else75 74 DO jj = 2, jpjm1 76 75 DO ji = 2, jpim1 77 # endif78 76 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels 79 77 ikbv = mbkv(ji,jj) … … 82 80 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 83 81 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 82 83 ! (ISF) stability criteria for top friction 84 ikbu = miku(ji,jj) ! first wet ocean u- & v-levels 85 ikbv = mikv(ji,jj) 86 ! 87 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 88 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( tfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) & 89 & * (1.-umask(ji,jj,1)) 90 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( tfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) & 91 & * (1.-vmask(ji,jj,1)) 92 ! (ISF) 93 84 94 END DO 85 95 END DO … … 89 99 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 90 100 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 91 CALL trd_ mod( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_trd_bfr, 'DYN', kt )101 CALL trd_dyn( ztrdu(:,:,:), ztrdv(:,:,:), jpdyn_bfr, kt ) 92 102 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 93 103 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4624 r4990 29 29 !!---------------------------------------------------------------------- 30 30 USE oce ! ocean dynamics and tracers 31 USE sbc_oce ! surface variable (only for the flag with ice shelf) 31 32 USE dom_oce ! ocean space and time domain 32 33 USE phycst ! physical constants 33 USE trdmod ! ocean dynamics trends 34 USE trdmod_oce ! ocean variables trends 34 USE trd_oce ! trends: ocean variables 35 USE trddyn ! trend manager: dynamics 36 ! 35 37 USE in_out_manager ! I/O manager 36 38 USE prtctl ! Print control 37 USE lbclnk ! lateral boundary condition 39 USE lbclnk ! lateral boundary condition 38 40 USE lib_mpp ! MPP library 41 USE eosbn2 ! compute density 39 42 USE wrk_nemo ! Memory Allocation 40 43 USE timing ! Timing … … 74 77 !! 75 78 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 76 !! - Save the trend(l_trddyn=T)79 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 77 80 !!---------------------------------------------------------------------- 78 81 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 99 102 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 100 103 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 101 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )104 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 102 105 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 103 106 ENDIF … … 175 178 IF( ln_hpg_prj ) ioptio = ioptio + 1 176 179 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 180 IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 ) & 181 & CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity.' ) 177 182 ! 178 183 END SUBROUTINE dyn_hpg_init … … 315 320 316 321 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) 317 # if defined key_vectopt_loop318 jj = 1319 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)320 # else321 322 DO jj = 2, jpjm1 322 323 DO ji = 2, jpim1 323 # endif324 324 iku = mbku(ji,jj) 325 325 ikv = mbkv(ji,jj) … … 338 338 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 339 339 ENDIF 340 # if ! defined key_vectopt_loop 341 END DO 342 # endif 340 END DO 343 341 END DO 344 342 ! … … 368 366 INTEGER, INTENT(in) :: kt ! ocean time-step index 369 367 !! 370 INTEGER :: ji, jj, jk ! dummy loop indices 371 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 372 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 373 !!---------------------------------------------------------------------- 374 ! 375 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 376 ! 377 IF( kt == nit000 ) THEN 368 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices 369 REAL(wp) :: zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1 ! temporary scalars 370 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj, zrhd 371 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 372 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 373 !!---------------------------------------------------------------------- 374 ! 375 CALL wrk_alloc( jpi,jpj, 2, ztstop) 376 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 377 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 378 ! 379 IF( kt == nit000 ) THEN 378 380 IF(lwp) WRITE(numout,*) 379 381 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' … … 384 386 zcoef0 = - grav * 0.5_wp 385 387 ! To use density and not density anomaly 386 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 387 ELSE ; znad = 0._wp ! Fixed volume 388 ENDIF 389 390 ! Surface value 388 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 389 ! ELSE ; znad = 0._wp ! Fixed volume 390 ! ENDIF 391 znad=1._wp 392 ! iniitialised to 0. zhpi zhpi 393 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 394 395 !================================================================================== 396 !=====Compute iceload and contribution of the half first wet layer ================= 397 !=================================================================================== 398 399 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 400 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 401 402 ! compute density of the water displaced by the ice shelf 403 zrhd = rhd ! save rhd 404 DO jk = 1, jpk 405 zdept(:,:)=gdept_1d(jk) 406 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk)) 407 END DO 408 WHERE ( tmask(:,:,:) == 1._wp) 409 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 410 END WHERE 411 412 ! compute rhd at the ice/oce interface (ice shelf side) 413 CALL eos(ztstop,risfdep,zrhdtop_isf) 414 415 ! compute rhd at the ice/oce interface (ocean side) 416 DO ji=1,jpi 417 DO jj=1,jpj 418 ikt=mikt(ji,jj) 419 ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 420 ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 421 END DO 422 END DO 423 CALL eos(ztstop,risfdep,zrhdtop_oce) 424 ! 425 ! Surface value + ice shelf gradient 426 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 427 ziceload = 0._wp 428 DO jj = 1, jpj 429 DO ji = 1, jpi ! vector opt. 430 ikt=mikt(ji,jj) 431 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 432 DO jk=2,ikt-1 433 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 434 & * (1._wp - tmask(ji,jj,jk)) 435 END DO 436 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 437 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 438 END DO 439 END DO 440 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 441 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 391 442 DO jj = 2, jpjm1 392 443 DO ji = fs_2, fs_jpim1 ! vector opt. 393 ! hydrostatic pressure gradient along s-surfaces 394 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 395 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 396 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 397 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 398 ! s-coordinate pressure gradient correction 444 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 445 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 446 ! we assume ISF is in isostatic equilibrium 447 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj ,iktp1i) & 448 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 449 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 450 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 451 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)) ) 452 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji ,jj+1,iktp1j) & 453 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 454 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 455 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 456 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) 457 ! s-coordinate pressure gradient correction (=0 if z coordinate) 399 458 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 400 459 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) … … 402 461 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 403 462 ! add to the general momentum trend 404 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 405 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 406 END DO 407 END DO 408 409 ! interior value (2=<jk=<jpkm1) 410 DO jk = 2, jpkm1 411 DO jj = 2, jpjm1 412 DO ji = fs_2, fs_jpim1 ! vector opt. 463 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 464 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 465 END DO 466 END DO 467 !================================================================================== 468 !===== Compute partial cell contribution for the top cell ========================= 469 !================================================================================== 470 DO jj = 2, jpjm1 471 DO ji = fs_2, fs_jpim1 ! vector opt. 472 iku = miku(ji,jj) ; 473 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 474 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 475 ! u direction 476 IF ( iku .GT. 1 ) THEN 477 ! case iku 478 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu & 479 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) & 480 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 481 ! corrective term ( = 0 if z coordinate ) 482 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 483 ! zhpi will be added in interior loop 484 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap 485 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 486 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku) 487 488 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 489 zhpiint = zcoef0 / e1u(ji,jj) & 490 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) & 491 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) & 492 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) & 493 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) ) 494 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 495 END IF 496 497 ! v direction 498 ikv = mikv(ji,jj) 499 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 500 IF ( ikv .GT. 1 ) THEN 501 ! case ikv 502 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv & 503 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) & 504 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 505 ! corrective term ( = 0 if z coordinate ) 506 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 507 ! zhpi will be added in interior loop 508 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap 509 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 510 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv) 511 512 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 513 zhpjint = zcoef0 / e2v(ji,jj) & 514 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) & 515 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) & 516 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) & 517 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) ) 518 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 519 END IF 520 END DO 521 END DO 522 523 !================================================================================== 524 !===== Compute interior value ===================================================== 525 !================================================================================== 526 527 DO jj = 2, jpjm1 528 DO ji = fs_2, fs_jpim1 ! vector opt. 529 iku=miku(ji,jj); ikv=mikv(ji,jj) 530 DO jk = 2, jpkm1 413 531 ! hydrostatic pressure gradient along s-surfaces 414 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 415 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 416 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 417 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 418 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 419 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 532 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 533 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 534 & + zcoef0 / e1u(ji,jj) & 535 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 536 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 537 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 538 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 420 539 ! s-coordinate pressure gradient correction 421 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 422 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 423 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 424 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 540 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 541 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 542 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 543 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 544 545 ! hydrostatic pressure gradient along s-surfaces 546 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 547 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 548 & + zcoef0 / e2v(ji,jj) & 549 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 550 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 551 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 552 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 553 ! s-coordinate pressure gradient correction 554 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 555 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 556 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 425 557 ! add to the general momentum trend 426 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 427 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 558 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 428 559 END DO 429 560 END DO 430 561 END DO 431 ! 432 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 562 563 !================================================================================== 564 !===== Compute bottom cell contribution (partial cell) ============================ 565 !================================================================================== 566 567 # if defined key_vectopt_loop 568 jj = 1 569 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 570 # else 571 DO jj = 2, jpjm1 572 DO ji = 2, jpim1 573 # endif 574 iku = mbku(ji,jj) 575 ikv = mbkv(ji,jj) 576 577 IF (iku .GT. 1) THEN 578 ! remove old value (interior case) 579 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 580 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 581 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 582 ! put new value 583 ! -zpshpi to avoid double contribution of the partial step in the top layer 584 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 585 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 586 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 587 END IF 588 ! v direction 589 IF (ikv .GT. 1) THEN 590 ! remove old value (interior case) 591 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 592 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 593 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 594 ! put new value 595 ! -zpshpj to avoid double contribution of the partial step in the top layer 596 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 597 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 598 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 599 END IF 600 # if ! defined key_vectopt_loop 601 END DO 602 # endif 603 END DO 604 605 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 606 rhd = zrhd 607 ! 608 CALL wrk_dealloc( jpi,jpj,2, ztstop) 609 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 610 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 433 611 ! 434 612 END SUBROUTINE hpg_sco 613 435 614 436 615 SUBROUTINE hpg_djc( kt ) … … 671 850 !! 672 851 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 673 !! - Save the trend (l_trddyn=T)674 !!675 852 !!---------------------------------------------------------------------- 676 853 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear … … 724 901 725 902 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 726 DO jj = 1, jpj; DO ji = 1, jpi 727 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 728 END DO ; END DO 729 730 DO jk = 2, jpk; DO jj = 1, jpj; DO ji = 1, jpi 731 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 732 END DO ; END DO ; END DO 733 734 fsp(:,:,:) = zrhh(:,:,:) 903 DO jj = 1, jpj 904 DO ji = 1, jpi 905 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 906 END DO 907 END DO 908 909 DO jk = 2, jpk 910 DO jj = 1, jpj 911 DO ji = 1, jpi 912 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 913 END DO 914 END DO 915 END DO 916 917 fsp(:,:,:) = zrhh (:,:,:) 735 918 xsp(:,:,:) = zdept(:,:,:) 736 919 … … 933 1116 END SUBROUTINE hpg_prj 934 1117 1118 935 1119 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 936 1120 !!---------------------------------------------------------------------- … … 940 1124 !! 941 1125 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 1126 !! 942 1127 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 943 !!944 1128 !!---------------------------------------------------------------------- 945 1129 IMPLICIT NONE … … 949 1133 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 950 1134 ! 2: Linear 951 952 ! Local Variables 1135 ! 953 1136 INTEGER :: ji, jj, jk ! dummy loop indices 954 1137 INTEGER :: jpi, jpj, jpkm1 … … 1040 1223 ENDIF 1041 1224 1042 1043 1225 END SUBROUTINE cspline 1044 1226 … … 1050 1232 !! ** Purpose : 1-d linear interpolation 1051 1233 !! 1052 !! ** Method : 1053 !! interpolation is straight forward 1234 !! ** Method : interpolation is straight forward 1054 1235 !! extrapolation is also permitted (no value limit) 1055 !!1056 1236 !!---------------------------------------------------------------------- 1057 1237 IMPLICIT NONE … … 1070 1250 END FUNCTION interp1 1071 1251 1252 1072 1253 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1073 1254 !!---------------------------------------------------------------------- … … 1133 1314 END FUNCTION integ_spline 1134 1315 1135 1136 1316 !!====================================================================== 1137 1317 END MODULE dynhpg -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90
r3294 r4990 14 14 USE oce ! ocean dynamics and tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE trdmod ! ocean dynamics trends 17 USE trdmod_oce ! ocean variables trends 16 USE trd_oce ! trends: ocean variables 17 USE trddyn ! trend manager: dynamics 18 ! 18 19 USE in_out_manager ! I/O manager 19 20 USE lib_mpp ! MPP library … … 52 53 !! 53 54 !! ** Action : - Update the (ua, va) with the hor. ke gradient trend 54 !! - s ave this trends(l_trddyn=T) for post-processing55 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 55 56 !!---------------------------------------------------------------------- 56 57 INTEGER, INTENT( in ) :: kt ! ocean time-step index 57 ! !58 ! 58 59 INTEGER :: ji, jj, jk ! dummy loop indices 59 60 REAL(wp) :: zu, zv ! temporary scalars … … 131 132 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 132 133 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 133 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_keg, 'DYN', kt )134 CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt ) 134 135 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 135 136 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf.F90
r4522 r4990 15 15 USE phycst ! physical constants 16 16 USE ldfdyn_oce ! ocean dynamics lateral physics 17 USE ldftra_oce ! ocean tracers lateral physics 17 18 USE ldfslp ! lateral mixing: slopes of mixing orientation 18 19 USE dynldf_bilapg ! lateral mixing (dyn_ldf_bilapg routine) … … 20 21 USE dynldf_iso ! lateral mixing (dyn_ldf_iso routine) 21 22 USE dynldf_lap ! lateral mixing (dyn_ldf_lap routine) 22 USE ldftra_oce, ONLY: ln_traldf_hor ! ocean tracers lateral physics23 USE trd mod ! ocean dynamics and tracer trends24 USE trdmod_oce ! ocean variables trends23 USE trd_oce ! trends: ocean variables 24 USE trddyn ! trend manager: dynamics (trd_dyn routine) 25 ! 25 26 USE prtctl ! Print control 26 27 USE in_out_manager ! I/O manager … … 30 31 USE timing ! Timing 31 32 32 33 33 IMPLICIT NONE 34 34 PRIVATE … … 55 55 !! ** Purpose : compute the lateral ocean dynamics physics. 56 56 !!---------------------------------------------------------------------- 57 !58 57 INTEGER, INTENT(in) :: kt ! ocean time-step index 59 58 ! … … 107 106 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 108 107 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 109 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_ldf, 'DYN', kt )108 CALL trd_dyn( ztrdu, ztrdv, jpdyn_ldf, kt ) 110 109 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 111 110 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90
r3634 r4990 19 19 USE dom_oce ! ocean space and time domain 20 20 USE ldfdyn_oce ! ocean dynamics: lateral physics 21 ! 21 22 USE in_out_manager ! I/O manager 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 25 24 USE wrk_nemo ! Memory Allocation … … 70 69 !! Add this before trend to the general trend (ua,va): 71 70 !! (ua,va) = (ua,va) + (diffu,diffv) 72 !! 'key_trddyn' defined: the two components of the horizontal73 !! diffusion trend are saved.74 71 !! 75 72 !! ** Action : - Update (ua,va) with the before iso-level biharmonic 76 73 !! mixing trend. 77 74 !!---------------------------------------------------------------------- 78 !79 75 INTEGER, INTENT(in) :: kt ! ocean time-step index 80 76 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90
r4488 r4990 19 19 USE dom_oce ! ocean space and time domain 20 20 USE ldfdyn_oce ! ocean dynamics lateral physics 21 USE zdf_oce ! ocean vertical physics 22 USE ldfslp ! iso-neutral slopes available 21 23 USE ldftra_oce, ONLY: ln_traldf_iso 22 USE zdf_oce ! ocean vertical physics 23 USE trdmod ! ocean dynamics trends 24 USE trdmod_oce ! ocean variables trends 25 USE ldfslp ! iso-neutral slopes available 24 ! 26 25 USE in_out_manager ! I/O manager 27 26 USE lib_mpp ! MPP library … … 81 80 !! -3- Add this trend to the general trend (ta,sa): 82 81 !! (ua,va) = (ua,va) + (zwk3,zwk4) 83 !! 'key_trddyn' defined: the trend is saved for diagnostics.84 82 !! 85 83 !! ** Action : - Update (ua,va) arrays with the before geopotential 86 84 !! biharmonic mixing trend. 87 !! - save the trend in (zwk3,zwk4) ('key_trddyn')88 85 !!---------------------------------------------------------------------- 89 86 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 201 198 !! pu and pv (all the components except 202 199 !! second order vertical derivative term) 203 !! 'key_trddyn' defined: the trend is saved for diagnostics. 204 !!---------------------------------------------------------------------- 205 !! 200 !!---------------------------------------------------------------------- 206 201 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity 207 202 ! ! 2nd call: ahm x these fields -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90
r4488 r4990 22 22 USE ldftra_oce ! ocean tracer lateral physics 23 23 USE zdf_oce ! ocean vertical physics 24 USE trdmod ! ocean dynamics trends25 USE trdmod_oce ! ocean variables trends26 24 USE ldfslp ! iso-neutral slopes 25 ! 27 26 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 28 27 USE in_out_manager ! I/O manager -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap.F90
r3294 r4990 19 19 USE ldfdyn_oce ! ocean dynamics: lateral physics 20 20 USE zdf_oce ! ocean vertical physics 21 ! 21 22 USE in_out_manager ! I/O manager 22 USE trdmod ! ocean dynamics trends23 USE trdmod_oce ! ocean variables trends24 USE ldfslp ! iso-neutral slopes25 23 USE timing ! Timing 26 24 … … 57 55 !! Add this before trend to the general trend (ua,va): 58 56 !! (ua,va) = (ua,va) + (diffu,diffv) 59 !! 'key_trddyn' activated: the two components of the horizontal60 !! diffusion trend are saved.61 57 !! 62 !! ** Action : - Update (ua,va) with the before iso-level harmonic 63 !! mixing trend. 58 !! ** Action : - Update (ua,va) with the iso-level harmonic mixing trend 64 59 !!---------------------------------------------------------------------- 65 60 INTEGER, INTENT( in ) :: kt ! ocean time-step index -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r4370 r4990 18 18 !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 19 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 20 !! 3.7 ! 2014-04 (G. Madec) add the diagnostic of the time filter trends 20 21 !!------------------------------------------------------------------------- 21 22 … … 34 35 USE bdydyn ! ocean open boundary conditions 35 36 USE bdyvol ! ocean open boundary condition (bdy_vol routines) 37 USE trd_oce ! trends: ocean variables 38 USE trddyn ! trend manager: dynamics 39 USE trdken ! trend manager: kinetic energy 40 ! 36 41 USE in_out_manager ! I/O manager 42 USE iom ! I/O manager library 37 43 USE lbclnk ! lateral boundary condition (or mpp link) 38 44 USE lib_mpp ! MPP library 39 45 USE wrk_nemo ! Memory Allocation 40 46 USE prtctl ! Print control 41 47 USE timing ! Timing 42 48 #if defined key_agrif 43 49 USE agrif_opa_interp 44 50 #endif 45 USE timing ! Timing46 51 47 52 IMPLICIT NONE … … 79 84 !! at the local domain boundaries through lbc_lnk call, 80 85 !! at the one-way open boundaries (lk_bdy=T), 81 !! at the AGRIF zoom 86 !! at the AGRIF zoom boundaries (lk_agrif=T) 82 87 !! 83 88 !! * Apply the time filter applied and swap of the dynamics … … 99 104 REAL(wp) :: z2dt ! temporary scalar 100 105 #endif 101 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars102 REAL(wp) :: zve3a, zve3n, zve3b, zvf 103 REAL(wp), POINTER, DIMENSION(:,:) :: zu a, zva104 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 107 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 108 REAL(wp), POINTER, DIMENSION(:,:) :: zue, zve 109 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva 105 110 !!---------------------------------------------------------------------- 106 111 ! 107 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt')108 ! 109 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f)110 IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva)112 IF( nn_timing == 1 ) CALL timing_start('dyn_nxt') 113 ! 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 115 IF( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zue, zve ) 111 116 ! 112 117 IF( kt == nit000 ) THEN … … 152 157 153 158 # if defined key_dynspg_ts 159 !!gm IF ( lk_dynspg_ts ) THEN .... 154 160 ! Ensure below that barotropic velocities match time splitting estimate 155 161 ! Compute actual transport and replace it with ts estimate at "after" time step 156 zu a(:,:) = 0._wp157 zv a(:,:) = 0._wp158 DO jk = 1, jpkm1159 zu a(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)160 zv a(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)162 zue(:,:) = fse3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1) 163 zve(:,:) = fse3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1) 164 DO jk = 2, jpkm1 165 zue(:,:) = zue(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zve(:,:) = zve(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 161 167 END DO 162 168 DO jk = 1, jpkm1 163 ua(:,:,jk) = ( ua(:,:,jk) - zu a(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)164 va(:,:,jk) = ( va(:,:,jk) - zv a(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)169 ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * hur_a(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * hvr_a(:,:) + va_b(:,:) ) * vmask(:,:,jk) 165 171 END DO 166 172 … … 175 181 END DO 176 182 ENDIF 183 !!gm ENDIF 177 184 # endif 178 185 … … 195 202 # endif 196 203 #endif 204 205 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 206 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step 207 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt 208 ! 209 ! ! Kinetic energy and Conversion 210 IF( ln_KE_trd ) CALL trd_dyn( ua, va, jpdyn_ken, kt ) 211 ! 212 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 213 zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt 214 zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt 215 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin time filter 216 CALL iom_put( "vtrd_tot", zva ) 217 ENDIF 218 ! 219 zua(:,:,:) = un(:,:,:) ! save the now velocity before the asselin filter 220 zva(:,:,:) = vn(:,:,:) ! (caution: there will be a shift by 1 timestep in the 221 ! ! computation of the asselin filter trends) 222 ENDIF 197 223 198 224 ! Time filter and swap of dynamics arrays … … 217 243 DO jj = 1, jpj 218 244 DO ji = 1, jpi 219 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2. e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) )220 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2. e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) )245 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 246 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 221 247 ! 222 248 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 301 327 ! Revert "before" velocities to time split estimate 302 328 ! Doing it here also means that asselin filter contribution is removed 303 zu a(:,:) = 0._wp304 zv a(:,:) = 0._wp305 DO jk = 1, jpkm1306 zu a(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)307 zv a(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)329 zue(:,:) = fse3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1) 330 zve(:,:) = fse3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1) 331 DO jk = 2, jpkm1 332 zue(:,:) = zue(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 333 zve(:,:) = zve(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 308 334 END DO 309 335 DO jk = 1, jpkm1 310 ub(:,:,jk) = ub(:,:,jk) - (zu a(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk)311 vb(:,:,jk) = vb(:,:,jk) - (zv a(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk)336 ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 337 vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 312 338 END DO 313 339 ENDIF … … 327 353 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 328 354 END DO 329 hur_b(:,:) = umask (:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )330 hvr_b(:,:) = vmask (:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )355 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 356 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 331 357 ENDIF 332 358 ! … … 335 361 ! 336 362 DO jk = 1, jpkm1 337 #if defined key_vectopt_loop338 DO jj = 1, 1 !Vector opt. => forced unrolling339 DO ji = 1, jpij340 #else341 363 DO jj = 1, jpj 342 364 DO ji = 1, jpi 343 #endif344 365 un_b(ji,jj) = un_b(ji,jj) + fse3u_a(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 345 366 vn_b(ji,jj) = vn_b(ji,jj) + fse3v_a(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) … … 358 379 ! 359 380 ! 381 382 IF( l_trddyn ) THEN ! 3D output: asselin filter trends on momentum 383 zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt 384 zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt 385 CALL trd_dyn( zua, zva, jpdyn_atf, kt ) 386 ENDIF 387 ! 360 388 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & 361 389 & tab3d_2=vn, clinfo2=' Vn: ' , mask2=vmask ) 362 390 ! 363 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f)364 IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva)391 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f, zua, zva ) 392 IF( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zue, zve ) 365 393 ! 366 394 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') … … 370 398 !!========================================================================= 371 399 END MODULE dynnxt 372 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r4496 r4990 26 26 USE sbctide 27 27 USE updtide 28 USE trdmod ! ocean dynamics trends 29 USE trdmod_oce ! ocean variables trends 28 USE trd_oce ! trends: ocean variables 29 USE trddyn ! trend manager: dynamics 30 ! 30 31 USE prtctl ! Print control (prt_ctl routine) 31 32 USE in_out_manager ! I/O manager 32 33 USE lib_mpp ! MPP library 33 USE solver 34 USE wrk_nemo 35 USE timing 34 USE solver ! solver initialization 35 USE wrk_nemo ! Memory Allocation 36 USE timing ! Timing 36 37 37 38 … … 163 164 END DO 164 165 END DO 165 END DO 166 END DO 167 168 !!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ???? 169 166 170 ENDIF 167 171 … … 191 195 CASE( 2 ) 192 196 z2dt = 2. * rdt 193 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt197 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 194 198 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / z2dt - ztrdu(:,:,:) 195 199 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / z2dt - ztrdv(:,:,:) 196 200 END SELECT 197 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_spg, 'DYN', kt )201 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt ) 198 202 ! 199 203 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) … … 246 250 IF( ( ioptio > 1 .AND. .NOT. lk_esopa ) .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 247 251 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme with a key cpp' ) 252 IF( ( lk_dynspg_ts .OR. lk_dynspg_exp ) .AND. nn_isf .NE. 0 ) & 253 & CALL ctl_stop( ' dynspg_ts and dynspg_exp not tested with ice shelf cavity ' ) 248 254 ! 249 255 IF( lk_esopa ) nspg = -1 -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r4328 r4990 19 19 USE sbc_oce ! surface boundary condition: ocean 20 20 USE phycst ! physical constants 21 ! 21 22 USE in_out_manager ! I/O manager 22 23 USE lib_mpp ! distributed memory computing library -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r4328 r4990 13 13 !! - ! 2006-08 (J.Chanut, A.Sellar) Calls to BDY routines. 14 14 !! 3.2 ! 2009-03 (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module 15 !! 3.7 ! 2014-04 (F. Roquet, G. Madec) add some trends diag 15 16 !!---------------------------------------------------------------------- 16 17 #if defined key_dynspg_flt || defined key_esopa … … 36 37 USE bdyvol ! ocean open boundary condition (bdy_vol routine) 37 38 USE cla ! cross land advection 39 USE trd_oce ! trends: ocean variables 40 USE trddyn ! trend manager: dynamics 41 ! 38 42 USE in_out_manager ! I/O manager 39 43 USE lib_mpp ! distributed memory computing library … … 43 47 USE iom 44 48 USE lib_fortran 49 USE timing ! Timing 45 50 #if defined key_agrif 46 51 USE agrif_opa_interp 47 52 #endif 48 USE timing ! Timing49 53 50 54 IMPLICIT NONE … … 99 103 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 100 104 !! 101 !! References : Roullet and Madec 1999, JGR.105 !! References : Roullet and Madec, JGR, 2000. 102 106 !!--------------------------------------------------------------------- 103 107 INTEGER, INTENT(in ) :: kt ! ocean time-step index 104 108 INTEGER, INTENT( out) :: kindic ! solver convergence flag (<0 if not converge) 105 ! !109 ! 106 110 INTEGER :: ji, jj, jk ! dummy loop indices 107 111 REAL(wp) :: z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv ! local scalars 112 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 113 REAL(wp), POINTER, DIMENSION(:,:) :: zpw 108 114 !!---------------------------------------------------------------------- 109 115 ! 110 116 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_flt') 111 !112 117 ! 113 118 IF( kt == nit000 ) THEN … … 179 184 END DO 180 185 ! 186 IF( l_trddyn ) THEN ! temporary save of spg trends 187 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 188 DO jk = 1, jpkm1 ! unweighted time stepping 189 DO jj = 2, jpjm1 190 DO ji = fs_2, fs_jpim1 ! vector opt. 191 ztrdu(ji,jj,jk) = spgu(ji,jj) * umask(ji,jj,jk) 192 ztrdv(ji,jj,jk) = spgv(ji,jj) * vmask(ji,jj,jk) 193 END DO 194 END DO 195 END DO 196 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgexp, kt ) 197 ENDIF 198 ! 181 199 ENDIF 182 200 … … 194 212 DO jj = 2, jpjm1 195 213 DO ji = fs_2, fs_jpim1 ! vector opt. 196 spgu(ji,jj) = 0._wp 197 spgv(ji,jj) = 0._wp 198 END DO 199 END DO 200 201 ! vertical sum 202 !CDIR NOLOOPCHG 203 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 204 DO jk = 1, jpkm1 205 DO ji = 1, jpij 206 spgu(ji,1) = spgu(ji,1) + fse3u_a(ji,1,jk) * ua(ji,1,jk) 207 spgv(ji,1) = spgv(ji,1) + fse3v_a(ji,1,jk) * va(ji,1,jk) 208 END DO 209 END DO 210 ELSE ! No vector opt. 211 DO jk = 1, jpkm1 212 DO jj = 2, jpjm1 213 DO ji = 2, jpim1 214 spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk) 215 spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk) 216 END DO 217 END DO 218 END DO 219 ENDIF 220 221 ! transport: multiplied by the horizontal scale factor 222 DO jj = 2, jpjm1 214 spgu(ji,jj) = fse3u_a(ji,jj,1) * ua(ji,jj,1) 215 spgv(ji,jj) = fse3v_a(ji,jj,1) * va(ji,jj,1) 216 END DO 217 END DO 218 DO jk = 2, jpkm1 ! vertical sum 219 DO jj = 2, jpjm1 220 DO ji = fs_2, fs_jpim1 ! vector opt. 221 spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk) 222 spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk) 223 END DO 224 END DO 225 END DO 226 227 DO jj = 2, jpjm1 ! transport: multiplied by the horizontal scale factor 223 228 DO ji = fs_2, fs_jpim1 ! vector opt. 224 229 spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj) … … 322 327 ENDIF 323 328 #endif 329 330 IF( l_trddyn ) THEN 331 ztrdu(:,:,:) = ua(:,:,:) ! save the after velocity before the filtered SPG 332 ztrdv(:,:,:) = va(:,:,:) 333 ! 334 CALL wrk_alloc( jpi, jpj, zpw ) 335 ! 336 zpw(:,:) = - z2dt * gcx(:,:) 337 CALL iom_put( "ssh_flt" , zpw ) ! output equivalent ssh modification due to implicit filter 338 ! 339 ! ! save surface pressure flux: -pw at z=0 340 zpw(:,:) = - rau0 * grav * sshn(:,:) * wn(:,:,1) * tmask(:,:,1) 341 CALL iom_put( "pw0_exp" , zpw ) 342 zpw(:,:) = wn(:,:,1) 343 CALL iom_put( "w0" , zpw ) 344 zpw(:,:) = rau0 * z2dtg * gcx(:,:) * wn(:,:,1) * tmask(:,:,1) 345 CALL iom_put( "pw0_flt" , zpw ) 346 ! 347 CALL wrk_dealloc( jpi, jpj, zpw ) 348 ! 349 ENDIF 350 324 351 ! Add the trends multiplied by z2dt to the after velocity 325 352 ! ------------------------------------------------------- … … 336 363 END DO 337 364 338 ! write filtered free surface arrays in restart file 339 ! -------------------------------------------------- 340 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) 341 ! 342 ! 343 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_flt') 365 IF( l_trddyn ) THEN ! save the explicit SPG trends for further diagnostics 366 ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / z2dt 367 ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / z2dt 368 CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgflt, kt ) 369 ! 370 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 371 ENDIF 372 373 IF( lrst_oce ) CALL flt_rst( kt, 'WRITE' ) ! write filtered free surface arrays in restart file 374 ! 375 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_flt') 344 376 ! 345 377 END SUBROUTINE dyn_spg_flt … … 352 384 !! ** Purpose : Read or write filtered free surface arrays in restart file 353 385 !!---------------------------------------------------------------------- 354 INTEGER , INTENT(in) :: kt 355 CHARACTER(len=*), INTENT(in) :: cdrw 386 INTEGER , INTENT(in) :: kt ! ocean time-step 387 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 356 388 !!---------------------------------------------------------------------- 357 389 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynvor.F90
r4624 r4990 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 18 !!---------------------------------------------------------------------- 18 19 … … 29 30 USE dommsk ! ocean mask 30 31 USE dynadv ! momentum advection (use ln_dynadv_vec value) 31 USE trd mod ! ocean dynamics trends32 USE trd mod_oce ! ocean variables trends32 USE trd_oce ! trends: ocean variables 33 USE trddyn ! trend manager: dynamics 33 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 34 35 USE prtctl ! Print control … … 73 74 !! ** Action : - Update (ua,va) with the now vorticity term trend 74 75 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative 75 !! and planetary vorticity trends) ('key_trddyn') 76 !! and planetary vorticity trends) and send them to trd_dyn 77 !! for futher diagnostics (l_trddyn=T) 76 78 !!---------------------------------------------------------------------- 77 79 INTEGER, INTENT( in ) :: kt ! ocean time-step index … … 108 110 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 109 111 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 110 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )112 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 111 113 ztrdu(:,:,:) = ua(:,:,:) 112 114 ztrdv(:,:,:) = va(:,:,:) … … 114 116 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 115 117 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 116 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 117 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 118 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 118 119 ELSE 119 120 CALL vor_ene( kt, ntot, ua, va ) ! total vorticity … … 127 128 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 128 129 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 129 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )130 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 130 131 ztrdu(:,:,:) = ua(:,:,:) 131 132 ztrdv(:,:,:) = va(:,:,:) … … 133 134 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 134 135 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 135 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 136 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 136 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 137 137 ELSE 138 138 CALL vor_ens( kt, ntot, ua, va ) ! total vorticity … … 146 146 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 147 147 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 148 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )148 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 149 149 ztrdu(:,:,:) = ua(:,:,:) 150 150 ztrdv(:,:,:) = va(:,:,:) … … 152 152 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 153 153 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 154 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 155 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 154 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 156 155 ELSE 157 156 CALL vor_mix( kt ) ! total vorticity (mix=ens-ene) … … 165 164 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 166 165 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 167 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_rvo, 'DYN', kt )166 CALL trd_dyn( ztrdu, ztrdv, jpdyn_rvo, kt ) 168 167 ztrdu(:,:,:) = ua(:,:,:) 169 168 ztrdv(:,:,:) = va(:,:,:) … … 171 170 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 172 171 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 173 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_pvo, 'DYN', kt ) 174 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_dat, 'DYN', kt ) 172 CALL trd_dyn( ztrdu, ztrdv, jpdyn_pvo, kt ) 175 173 ELSE 176 174 CALL vor_een( kt, ntot, ua, va ) ! total vorticity … … 211 209 !! 212 210 !! ** Action : - Update (ua,va) with the now vorticity term trend 213 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative214 !! and planetary vorticity trends) ('key_trddyn')215 211 !! 216 212 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 328 324 !! 329 325 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 330 !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative331 !! and planetary vorticity trends) ('key_trddyn')332 326 !! 333 327 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 444 438 !! 445 439 !! ** Action : - Update (ua,va) arrays with the now vorticity term trend 446 !! - Save the trends in (ztrdu,ztrdv) in 2 parts (relative447 !! and planetary vorticity trends) ('key_trddyn')448 440 !! 449 441 !! References : Sadourny, r., 1975, j. atmos. sciences, 32, 680-689. … … 557 549 !! 558 550 !! ** Action : - Update (ua,va) with the now vorticity term trend 559 !! - save the trends in (ztrdu,ztrdv) in 2 parts (relative560 !! and planetary vorticity trends) ('key_trddyn')561 551 !! 562 552 !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36 … … 601 591 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'dyn:vor_een : unable to allocate arrays' ) 602 592 ENDIF 603 ze3f(:,:,:) = 0. d0593 ze3f(:,:,:) = 0._wp 604 594 #endif 605 595 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r3294 r4990 16 16 USE dom_oce ! ocean space and time domain 17 17 USE sbc_oce ! surface boundary condition: ocean 18 USE trdmod_oce ! ocean variables trends 19 USE trdmod ! ocean dynamics trends 18 USE trd_oce ! trends: ocean variables 19 USE trddyn ! trend manager: dynamics 20 ! 20 21 USE in_out_manager ! I/O manager 21 USE lib_mpp 22 USE lib_mpp ! MPP library 22 23 USE prtctl ! Print control 23 USE wrk_nemo 24 USE timing 24 USE wrk_nemo ! Memory Allocation 25 USE timing ! Timing 25 26 26 27 IMPLICIT NONE 27 28 PRIVATE 28 29 29 PUBLIC dyn_zad ! routine called by step.F90 30 PUBLIC dyn_zad ! routine called by dynadv.F90 31 PUBLIC dyn_zad_zts ! routine called by dynadv.F90 30 32 31 33 !! * Substitutions … … 53 55 !! 54 56 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 55 !! - S ave the trends in (ztrdu,ztrdv) ('key_trddyn')57 !! - Send the trends to trddyn for diagnostics (l_trddyn=T) 56 58 !!---------------------------------------------------------------------- 57 59 INTEGER, INTENT(in) :: kt ! ocean time-step inedx … … 83 85 DO jj = 2, jpj ! vertical fluxes 84 86 DO ji = fs_2, jpi ! vector opt. 85 zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)87 zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 86 88 END DO 87 89 END DO … … 95 97 DO jj = 2, jpjm1 ! Surface and bottom values set to zero 96 98 DO ji = fs_2, fs_jpim1 ! vector opt. 97 zwuw(ji,jj, 1 ) = 0.e098 zwvw(ji,jj, 1 ) = 0.e099 zwuw(ji,jj,jpk) = 0. e0100 zwvw(ji,jj,jpk) = 0. e099 zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 100 zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 101 zwuw(ji,jj,jpk) = 0._wp 102 zwvw(ji,jj,jpk) = 0._wp 101 103 END DO 102 104 END DO … … 118 120 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 119 121 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 120 CALL trd_ mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt)122 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 121 123 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 122 124 ENDIF … … 132 134 END SUBROUTINE dyn_zad 133 135 136 SUBROUTINE dyn_zad_zts ( kt ) 137 !!---------------------------------------------------------------------- 138 !! *** ROUTINE dynzad_zts *** 139 !! 140 !! ** Purpose : Compute the now vertical momentum advection trend and 141 !! add it to the general trend of momentum equation. This version 142 !! uses sub-timesteps for improved numerical stability with small 143 !! vertical grid sizes. This is especially relevant when using 144 !! embedded ice with thin surface boxes. 145 !! 146 !! ** Method : The now vertical advection of momentum is given by: 147 !! w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] 148 !! w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ] 149 !! Add this trend to the general trend (ua,va): 150 !! (ua,va) = (ua,va) + w dz(u,v) 151 !! 152 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 153 !! - Save the trends in (ztrdu,ztrdv) ('key_trddyn') 154 !!---------------------------------------------------------------------- 155 INTEGER, INTENT(in) :: kt ! ocean time-step inedx 156 ! 157 INTEGER :: ji, jj, jk, jl ! dummy loop indices 158 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 159 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 160 REAL(wp) :: zua, zva ! temporary scalars 161 REAL(wp) :: zr_rdt ! temporary scalar 162 REAL(wp) :: z2dtzts ! length of Euler forward sub-timestep for vertical advection 163 REAL(wp) :: zts ! length of sub-timestep for vertical advection 164 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw, zww 165 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 166 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zus , zvs 167 !!---------------------------------------------------------------------- 168 ! 169 IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts') 170 ! 171 CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 172 CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 173 ! 174 IF( kt == nit000 ) THEN 175 IF(lwp)WRITE(numout,*) 176 IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps' 177 ENDIF 178 179 IF( l_trddyn ) THEN ! Save ua and va trends 180 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 181 ztrdu(:,:,:) = ua(:,:,:) 182 ztrdv(:,:,:) = va(:,:,:) 183 ENDIF 184 185 IF( neuler == 0 .AND. kt == nit000 ) THEN 186 z2dtzts = rdt / REAL( jnzts, wp ) ! = rdt (restart with Euler time stepping) 187 ELSE 188 z2dtzts = 2._wp * rdt / REAL( jnzts, wp ) ! = 2 rdt (leapfrog) 189 ENDIF 190 191 DO jk = 2, jpkm1 ! Calculate and store vertical fluxes 192 DO jj = 2, jpj 193 DO ji = fs_2, jpi ! vector opt. 194 zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 195 END DO 196 END DO 197 END DO 198 199 DO jj = 2, jpjm1 ! Surface and bottom advective fluxes set to zero 200 DO ji = fs_2, fs_jpim1 ! vector opt. 201 zwuw(ji,jj, 1:miku(ji,jj) ) = 0._wp 202 zwvw(ji,jj, 1:mikv(ji,jj) ) = 0._wp 203 zwuw(ji,jj,jpk) = 0._wp 204 zwvw(ji,jj,jpk) = 0._wp 205 END DO 206 END DO 207 208 ! Start with before values and use sub timestepping to reach after values 209 210 zus(:,:,:,1) = ub(:,:,:) 211 zvs(:,:,:,1) = vb(:,:,:) 212 213 DO jl = 1, jnzts ! Start of sub timestepping loop 214 215 IF( jl == 1 ) THEN ! Euler forward to kick things off 216 jtb = 1 ; jtn = 1 ; jta = 2 217 zts = z2dtzts 218 ELSEIF( jl == 2 ) THEN ! First leapfrog step 219 jtb = 1 ; jtn = 2 ; jta = 3 220 zts = 2._wp * z2dtzts 221 ELSE ! Shuffle pointers for subsequent leapfrog steps 222 jtb = MOD(jtb,3) + 1 223 jtn = MOD(jtn,3) + 1 224 jta = MOD(jta,3) + 1 225 ENDIF 226 227 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 228 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) 231 zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) 232 END DO 233 END DO 234 END DO 235 DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points 236 DO jj = 2, jpjm1 237 DO ji = fs_2, fs_jpim1 ! vector opt. 238 ! ! vertical momentum advective trends 239 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 240 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 241 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 242 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts 243 END DO 244 END DO 245 END DO 246 247 END DO ! End of sub timestepping loop 248 249 zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts ) 250 DO jk = 1, jpkm1 ! Recover trends over the outer timestep 251 DO jj = 2, jpjm1 252 DO ji = fs_2, fs_jpim1 ! vector opt. 253 ! ! vertical momentum advective trends 254 ! ! add the trends to the general momentum trends 255 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt 256 va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt 257 END DO 258 END DO 259 END DO 260 261 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic 262 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 263 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 264 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zad, kt ) 265 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 266 ENDIF 267 ! ! Control print 268 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, & 269 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 270 ! 271 CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 272 CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 273 ! 274 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts') 275 ! 276 END SUBROUTINE dyn_zad_zts 277 134 278 !!====================================================================== 135 279 END MODULE dynzad -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r3294 r4990 20 20 21 21 USE ldfdyn_oce ! ocean dynamics: lateral physics 22 USE trd mod ! ocean active dynamics and tracers trends23 USE trd mod_oce ! ocean variables trends22 USE trd_oce ! trends: ocean variables 23 USE trddyn ! trend manager: dynamics 24 24 USE in_out_manager ! I/O manager 25 25 USE lib_mpp ! MPP library … … 91 91 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 92 92 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 93 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_zdf, 'DYN', kt ) 94 ! 93 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 95 94 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 96 95 ENDIF -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r4370 r4990 70 70 REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars 71 71 REAL(wp) :: ze3ua, ze3va 72 !!----------------------------------------------------------------------73 74 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws 75 73 !!---------------------------------------------------------------------- … … 101 99 102 100 IF( ln_bfrimp ) THEN 103 # if defined key_vectopt_loop104 DO jj = 1, 1105 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)106 # else107 101 DO jj = 2, jpjm1 108 102 DO ji = 2, jpim1 109 # endif110 103 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 111 104 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 112 105 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 113 106 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 107 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 108 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 109 IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu) 110 IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv) 114 111 END DO 115 112 END DO … … 138 135 ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 139 136 va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 140 END DO141 ! Add bottom stress due to barotropic component only:137 END DO 138 ! Add bottom/top stress due to barotropic component only: 142 139 DO jj = 2, jpjm1 143 140 DO ji = fs_2, fs_jpim1 ! vector opt. … … 148 145 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 149 146 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 147 ikbu = miku(ji,jj) ! top ocean level at u- and v-points 148 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 149 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu) 150 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv) 151 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 152 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va 150 153 END DO 151 154 END DO … … 166 169 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 167 170 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 168 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 171 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 169 172 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 170 173 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws … … 194 197 !----------------------------------------------------------------------- 195 198 ! 196 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 197 DO jj = 2, jpjm1 198 DO ji = fs_2, fs_jpim1 ! vector opt. 199 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 200 DO jj = 2, jpjm1 201 DO ji = fs_2, fs_jpim1 ! vector opt. 202 DO jk = miku(ji,jj)+1, jpkm1 199 203 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 200 204 END DO … … 204 208 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 205 209 DO ji = fs_2, fs_jpim1 ! vector opt. 206 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj, 1) + r_vvl * fse3u_a(ji,jj,1)210 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,miku(ji,jj)) + r_vvl * fse3u_a(ji,jj,miku(ji,jj)) 207 211 #if defined key_dynspg_ts 208 ua(ji,jj, 1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &212 ua(ji,jj,miku(ji,jj)) = ua(ji,jj,miku(ji,jj)) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 209 213 & / ( ze3ua * rau0 ) 210 214 #else 211 ua(ji,jj,1) = ub(ji,jj,1) + p2dt *(ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 212 & / ( fse3u(ji,jj,1) * rau0 ) ) 213 #endif 214 END DO 215 END DO 216 DO jk = 2, jpkm1 217 DO jj = 2, jpjm1 218 DO ji = fs_2, fs_jpim1 ! vector opt. 215 ua(ji,jj,miku(ji,jj)) = ub(ji,jj,miku(ji,jj)) & 216 & + p2dt *(ua(ji,jj,miku(ji,jj)) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 217 & / ( fse3u(ji,jj,miku(ji,jj)) * rau0 ) ) 218 #endif 219 DO jk = miku(ji,jj)+1, jpkm1 219 220 #if defined key_dynspg_ts 220 221 zrhs = ua(ji,jj,jk) ! zrhs=right hand side … … 230 231 DO ji = fs_2, fs_jpim1 ! vector opt. 231 232 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 232 END DO 233 END DO 234 DO jk = jpk-2, 1, -1 235 DO jj = 2, jpjm1 236 DO ji = fs_2, fs_jpim1 ! vector opt. 233 DO jk = jpk-2, miku(ji,jj), -1 237 234 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 238 235 END DO … … 292 289 !----------------------------------------------------------------------- 293 290 ! 294 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 295 DO jj = 2, jpjm1 296 DO ji = fs_2, fs_jpim1 ! vector opt. 291 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 292 DO jj = 2, jpjm1 293 DO ji = fs_2, fs_jpim1 ! vector opt. 294 DO jk = mikv(ji,jj)+1, jpkm1 297 295 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 298 296 END DO … … 302 300 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 303 301 DO ji = fs_2, fs_jpim1 ! vector opt. 304 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj, 1) + r_vvl * fse3v_a(ji,jj,1)302 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,mikv(ji,jj)) + r_vvl * fse3v_a(ji,jj,mikv(ji,jj)) 305 303 #if defined key_dynspg_ts 306 va(ji,jj, 1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &304 va(ji,jj,mikv(ji,jj)) = va(ji,jj,mikv(ji,jj)) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 307 305 & / ( ze3va * rau0 ) 308 306 #else 309 va(ji,jj,1) = vb(ji,jj,1) + p2dt *(va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 310 & / ( fse3v(ji,jj,1) * rau0 ) ) 311 #endif 312 END DO 313 END DO 314 DO jk = 2, jpkm1 315 DO jj = 2, jpjm1 316 DO ji = fs_2, fs_jpim1 ! vector opt. 307 va(ji,jj,mikv(ji,jj)) = vb(ji,jj,mikv(ji,jj)) & 308 & + p2dt *(va(ji,jj,mikv(ji,jj)) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 309 & / ( fse3v(ji,jj,mikv(ji,jj)) * rau0 ) ) 310 #endif 311 DO jk = mikv(ji,jj)+1, jpkm1 317 312 #if defined key_dynspg_ts 318 313 zrhs = va(ji,jj,jk) ! zrhs=right hand side … … 328 323 DO ji = fs_2, fs_jpim1 ! vector opt. 329 324 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 330 END DO 331 END DO 332 DO jk = jpk-2, 1, -1 333 DO jj = 2, jpjm1 334 DO ji = fs_2, fs_jpim1 ! vector opt. 325 DO jk = jpk-2, mikv(ji,jj), -1 335 326 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 336 327 END DO … … 352 343 !! restore bottom layer avmu(v) 353 344 IF( ln_bfrimp ) THEN 354 # if defined key_vectopt_loop 355 DO jj = 1, 1 356 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 357 # else 358 DO jj = 2, jpjm1 359 DO ji = 2, jpim1 360 # endif 361 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 362 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 363 avmu(ji,jj,ikbu+1) = 0.e0 364 avmv(ji,jj,ikbv+1) = 0.e0 365 END DO 366 END DO 345 DO jj = 2, jpjm1 346 DO ji = 2, jpim1 347 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 348 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 349 avmu(ji,jj,ikbu+1) = 0.e0 350 avmv(ji,jj,ikbv+1) = 0.e0 351 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 352 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 353 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0 354 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0 355 END DO 356 END DO 367 357 ENDIF 368 358 ! -
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r4486 r4990 31 31 USE bdy_par 32 32 USE bdydyn2d ! bdy_ssh routine 33 USE diaar5, ONLY: lk_diaar534 33 USE iom 35 34 #if defined key_agrif … … 111 110 ! 112 111 z1_rau0 = 0.5_wp * r1_rau0 113 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1)112 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 114 113 115 114 #if ! defined key_dynspg_ts … … 138 137 ! ! outputs ! 139 138 ! !------------------------------! 140 CALL iom_put( "ssh" , sshn 141 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height139 CALL iom_put( "ssh" , sshn ) ! sea surface height 140 if( iom_use('ssh2') ) CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 142 141 ! 143 142 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) … … 233 232 ! !------------------------------! 234 233 CALL iom_put( "woce", wn ) ! vertical velocity 235 IF( lk_diaar5 ) THEN! vertical mass transport & its square value234 IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value 236 235 CALL wrk_alloc( jpi, jpj, z2d ) 237 236 CALL wrk_alloc( jpi, jpj, jpk, z3d ) … … 241 240 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 242 241 END DO 243 CALL iom_put( "w_masstr" , z3d 244 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )242 CALL iom_put( "w_masstr" , z3d ) 243 IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 245 244 CALL wrk_dealloc( jpi, jpj, z2d ) 246 245 CALL wrk_dealloc( jpi, jpj, jpk, z3d ) … … 291 290 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 292 291 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 293 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)292 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * ssmask(:,:) 294 293 sshn(:,:) = ssha(:,:) ! now <-- after 295 294 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
