Changeset 12489 for NEMO/trunk/src/OCE/DYN
- Timestamp:
- 2020-02-28T16:55:11+01:00 (4 years ago)
- Location:
- NEMO/trunk/src/OCE/DYN
- Files:
-
- 7 edited
-
dynatf.F90 (modified) (9 diffs)
-
dynspg.F90 (modified) (6 diffs)
-
dynspg_exp.F90 (modified) (1 diff)
-
dynspg_ts.F90 (modified) (29 diffs)
-
dynzdf.F90 (modified) (13 diffs)
-
sshwzv.F90 (modified) (13 diffs)
-
wet_dry.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/DYN/dynatf.F90
r12377 r12489 87 87 !! arrays to start the next time step: 88 88 !! (puu(Kmm),pvv(Kmm)) = (puu(Kmm),pvv(Kmm)) 89 !! + atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ]89 !! + rn_atfp [ (puu(Kbb),pvv(Kbb)) + (puu(Kaa),pvv(Kaa)) - 2 (puu(Kmm),pvv(Kmm)) ] 90 90 !! Note that with flux form advection and non linear free surface, 91 91 !! the time filter is applied on thickness weighted velocity. … … 157 157 ! 158 158 IF( l_trddyn ) THEN ! prepare the atf trend computation + some diagnostics 159 z1_2dt = 1._wp / (2. * rdt) ! Euler or leap-frog time step160 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1._wp / rdt161 159 ! 162 160 ! ! Kinetic energy and Conversion … … 164 162 ! 165 163 IF( ln_dyn_trd ) THEN ! 3D output: total momentum trends 166 zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * z1_2dt167 zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * z1_2dt164 zua(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) * r1_Dt 165 zva(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) * r1_Dt 168 166 CALL iom_put( "utrd_tot", zua ) ! total momentum trends, except the asselin time filter 169 167 CALL iom_put( "vtrd_tot", zva ) … … 178 176 ! ------------------------------------------ 179 177 180 IF( .NOT. ( neuler == 0 .AND. kt == nit000 )) THEN !* Leap-Frog : Asselin time filter178 IF( .NOT. l_1st_euler ) THEN !* Leap-Frog : Asselin time filter 181 179 ! ! =============! 182 180 IF( ln_linssh ) THEN ! Fixed volume ! 183 181 ! ! =============! 184 182 DO_3D_11_11( 1, jpkm1 ) 185 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )186 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )183 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 184 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) 187 185 END_3D 188 186 ! ! ================! … … 193 191 ALLOCATE( ze3t_f(jpi,jpj,jpk), zwfld(jpi,jpj) ) 194 192 DO jk = 1, jpkm1 195 ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) )193 ze3t_f(:,:,jk) = pe3t(:,:,jk,Kmm) + rn_atfp * ( pe3t(:,:,jk,Kbb) - 2._wp * pe3t(:,:,jk,Kmm) + pe3t(:,:,jk,Kaa) ) 196 194 END DO 197 195 ! Add volume filter correction: compatibility with tracer advection scheme 198 196 ! => time filter + conservation correction 199 zcoef = atfp * rdt * r1_rau0197 zcoef = rn_atfp * rn_Dt * r1_rho0 200 198 zwfld(:,:) = emp_b(:,:) - emp(:,:) 201 199 IF ( ln_rnf ) zwfld(:,:) = zwfld(:,:) - ( rnf_b(:,:) - rnf(:,:) ) … … 209 207 ! to manage rnf, isf and possibly in the futur icb, tide water glacier (...) 210 208 ! ...(kt, coef, ktop, kbot, hz, fwf_b, fwf) 211 IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, atfp * rdt )209 IF ( ln_isf ) CALL isf_dynatf( kt, Kmm, ze3t_f, rn_atfp * rn_Dt ) 212 210 ! 213 211 pe3t(:,:,1:jpkm1,Kmm) = ze3t_f(:,:,1:jpkm1) ! filtered scale factor at T-points … … 218 216 CALL dom_vvl_interpol( pe3t(:,:,:,Kmm), pe3v(:,:,:,Kmm), 'V' ) 219 217 DO_3D_11_11( 1, jpkm1 ) 220 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) )221 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) )218 puu(ji,jj,jk,Kmm) = puu(ji,jj,jk,Kmm) + rn_atfp * ( puu(ji,jj,jk,Kbb) - 2._wp * puu(ji,jj,jk,Kmm) + puu(ji,jj,jk,Kaa) ) 219 pvv(ji,jj,jk,Kmm) = pvv(ji,jj,jk,Kmm) + rn_atfp * ( pvv(ji,jj,jk,Kbb) - 2._wp * pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk,Kaa) ) 222 220 END_3D 223 221 ! … … 236 234 zve3b = pe3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb) 237 235 ! 238 puu(ji,jj,jk,Kmm) = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk)239 pvv(ji,jj,jk,Kmm) = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk)236 puu(ji,jj,jk,Kmm) = ( zue3n + rn_atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 237 pvv(ji,jj,jk,Kmm) = ( zve3n + rn_atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 240 238 END_3D 241 239 pe3u(:,:,1:jpkm1,Kmm) = ze3u_f(:,:,1:jpkm1) … … 263 261 ENDIF 264 262 ! 265 ENDIF ! neuler /= 0263 ENDIF ! .NOT. l_1st_euler 266 264 ! 267 265 ! Set "now" and "before" barotropic velocities for next time step: -
NEMO/trunk/src/OCE/DYN/dynspg.F90
r12377 r12489 67 67 !! ln_apr_dyn=T : the atmospheric pressure forcing is applied 68 68 !! as the gradient of the inverse barometer ssh: 69 !! apgu = - 1/r au0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]70 !! apgv = - 1/r au0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]71 !! Note that as all external forcing a time averaging over a two r dt69 !! apgu = - 1/rho0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb] 70 !! apgv = - 1/rho0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb] 71 !! Note that as all external forcing a time averaging over a two rn_Dt 72 72 !! period is used to prevent the divergence of odd and even time step. 73 73 !!---------------------------------------------------------------------- … … 78 78 ! 79 79 INTEGER :: ji, jj, jk ! dummy loop indices 80 REAL(wp) :: z2dt, zg_2, zintp, zgr au0r, zld ! local scalars80 REAL(wp) :: z2dt, zg_2, zintp, zgrho0r, zld ! local scalars 81 81 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice 82 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv … … 114 114 ! 115 115 ! Update tide potential at the beginning of current time step 116 zt0step = REAL(nsec_day, wp)-0.5_wp*r dt116 zt0step = REAL(nsec_day, wp)-0.5_wp*rn_Dt 117 117 CALL upd_tide(zt0step, Kmm) 118 118 ! … … 134 134 ALLOCATE( zpice(jpi,jpj) ) 135 135 zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc ) 136 zgr au0r = - grav * r1_rau0137 zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zgr au0r136 zgrho0r = - grav * r1_rho0 137 zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zgrho0r 138 138 DO_2D_00_00 139 139 spgu(ji,jj) = spgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) … … 183 183 NAMELIST/namdyn_spg/ ln_dynspg_exp , ln_dynspg_ts, & 184 184 & ln_bt_fw, ln_bt_av , ln_bt_auto , & 185 & nn_ baro, rn_bt_cmax, nn_bt_flt, rn_bt_alpha185 & nn_e , rn_bt_cmax, nn_bt_flt, rn_bt_alpha 186 186 !!---------------------------------------------------------------------- 187 187 ! … … 222 222 ! 223 223 IF( nspg == np_TS ) THEN ! split-explicit scheme initialisation 224 CALL dyn_spg_ts_init ! do it first: set nn_ baroused to allocate some arrays later on224 CALL dyn_spg_ts_init ! do it first: set nn_e used to allocate some arrays later on 225 225 ENDIF 226 226 ! -
NEMO/trunk/src/OCE/DYN/dynspg_exp.F90
r12377 r12489 49 49 !! momentum trend the surface pressure gradient : 50 50 !! (uu(rhs),vv(rhs)) = (uu(rhs),vv(rhs)) + (spgu,spgv) 51 !! where spgu = -1/r au0 d/dx(ps) = -g/e1u di( ssh(now) )52 !! spgv = -1/r au0 d/dy(ps) = -g/e2v dj( ssh(now) )51 !! where spgu = -1/rho0 d/dx(ps) = -g/e1u di( ssh(now) ) 52 !! spgv = -1/rho0 d/dy(ps) = -g/e2v dj( ssh(now) ) 53 53 !! 54 54 !! ** Action : (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) trend of horizontal velocity increased by -
NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
r12377 r12489 72 72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 73 73 ! 74 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_ baro <= 2.5 nn_baro75 REAL(wp),SAVE :: r dtbt! Barotropic time step74 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e 75 REAL(wp),SAVE :: rDt_e ! Barotropic time step 76 76 ! 77 77 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields … … 102 102 ierr(:) = 0 103 103 ! 104 ALLOCATE( wgtbtp1(3*nn_ baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )104 ALLOCATE( wgtbtp1(3*nn_e), wgtbtp2(3*nn_e), zwz(jpi,jpj), STAT=ierr(1) ) 105 105 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 106 106 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) … … 150 150 LOGICAL :: ll_init ! =T : special startup of 2d equations 151 151 INTEGER :: noffset ! local integers : time offset for bdy update 152 REAL(wp) :: r1_ 2dt_b, z1_hu, z1_hv ! local scalars152 REAL(wp) :: r1_Dt_b, z1_hu, z1_hv ! local scalars 153 153 REAL(wp) :: za0, za1, za2, za3 ! - - 154 154 REAL(wp) :: zztmp, zldg ! - - … … 180 180 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 181 181 ! ! inverse of baroclinic time step 182 IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt ) 183 ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt ) 184 ENDIF 182 r1_Dt_b = 1._wp / rDt 185 183 ! 186 184 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 185 ll_fw_start = .FALSE. 188 186 ! ! time offset in steps for bdy data update 189 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_ baro187 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e 190 188 ELSE ; noffset = 0 191 189 ENDIF … … 198 196 IF(lwp) WRITE(numout,*) 199 197 ! 200 IF( neuler == 0) ll_init=.TRUE.201 ! 202 IF( ln_bt_fw .OR. neuler == 0) THEN198 IF( l_1st_euler ) ll_init=.TRUE. 199 ! 200 IF( ln_bt_fw .OR. l_1st_euler ) THEN 203 201 ll_fw_start =.TRUE. 204 202 noffset = 0 … … 209 207 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 210 208 ! 211 ENDIF 212 ! 213 ! If forward start at previous time step, and centered integration, 214 ! then update averaging weights: 215 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 216 ll_fw_start=.FALSE. 217 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 218 ENDIF 219 ! 220 209 ELSEIF( kt == nit000 + 1 ) THEN !* initialisation 2nd time-step 210 ! 211 IF( .NOT.ln_bt_fw ) THEN 212 ! If we did an Euler timestep on the first timestep we need to reset ll_fw_start 213 ! and the averaging weights. We don't have an easy way of telling whether we did 214 ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false. 215 ! at the end of the first timestep) so just do this in all cases. 216 ll_fw_start = .FALSE. 217 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 218 ENDIF 219 ! 220 ENDIF 221 ! 221 222 ! ----------------------------------------------------------------------------- 222 223 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) … … 302 303 IF( ln_bt_fw ) THEN ! Add wind forcing 303 304 DO_2D_00_00 304 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_r au0 * utau(ji,jj) * r1_hu(ji,jj,Kmm)305 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_r au0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm)305 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 306 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) 306 307 END_2D 307 308 ELSE 308 zztmp = r1_r au0 * r1_2309 zztmp = r1_rho0 * r1_2 309 310 DO_2D_00_00 310 311 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) … … 319 320 ! ! --------------------------------------------------- ! 320 321 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 321 zssh_frc(:,:) = r1_r au0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )322 zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 322 323 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 323 zztmp = r1_r au0 * r1_2324 zztmp = r1_rho0 * r1_2 324 325 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) & 325 326 & - rnf(:,:) - rnf_b(:,:) & … … 428 429 ! Update tide potential at the beginning of current time substep 429 430 IF( ln_tide_pot .AND. ln_tide ) THEN 430 zt0substep = REAL(nsec_day, wp) - 0.5_wp*r dt + (jn + noffset - 1) * rdt / REAL(nn_baro, wp)431 zt0substep = REAL(nsec_day, wp) - 0.5_wp*rn_Dt + (jn + noffset - 1) * rn_Dt / REAL(nn_e, wp) 431 432 CALL upd_tide(zt0substep, Kmm) 432 433 END IF … … 494 495 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV ) 495 496 #endif 496 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, r dtbt) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV497 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV 497 498 498 499 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where … … 509 510 DO_2D_00_00 510 511 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 511 ssha_e(ji,jj) = ( sshn_e(ji,jj) - r dtbt* ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj)512 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 512 513 END_2D 513 514 ! … … 599 600 DO_2D_00_00 600 601 ua_e(ji,jj) = ( un_e(ji,jj) & 601 & + r dtbt* ( zu_spg(ji,jj) &602 & + rDt_e * ( zu_spg(ji,jj) & 602 603 & + zu_trd(ji,jj) & 603 604 & + zu_frc(ji,jj) ) & … … 605 606 606 607 va_e(ji,jj) = ( vn_e(ji,jj) & 607 & + r dtbt* ( zv_spg(ji,jj) &608 & + rDt_e * ( zv_spg(ji,jj) & 608 609 & + zv_trd(ji,jj) & 609 610 & + zv_frc(ji,jj) ) & … … 624 625 ! 625 626 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 626 & + r dtbt* ( zhu_bck * zu_spg (ji,jj) & !627 & + rDt_e * ( zhu_bck * zu_spg (ji,jj) & ! 627 628 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 628 629 & + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu 629 630 ! 630 631 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 631 & + r dtbt* ( zhv_bck * zv_spg (ji,jj) & !632 & + rDt_e * ( zhv_bck * zv_spg (ji,jj) & ! 632 633 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 633 634 & + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv … … 637 638 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 638 639 DO_2D_00_00 639 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - r dtbt* zCdU_u(ji,jj) * hur_e(ji,jj))640 va_e(ji,jj) = va_e(ji,jj) /(1.0 - r dtbt* zCdU_v(ji,jj) * hvr_e(ji,jj))640 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 641 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) 641 642 END_2D 642 643 ENDIF … … 701 702 ! Set advection velocity correction: 702 703 IF (ln_bt_fw) THEN 703 IF( .NOT.( kt == nit000 .AND. neuler==0) ) THEN704 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 704 705 DO_2D_11_11 705 706 zun_save = un_adv(ji,jj) 706 707 zvn_save = vn_adv(ji,jj) 707 708 ! ! apply the previously computed correction 708 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) )709 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) )709 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) ) 710 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) ) 710 711 ! ! Update corrective fluxes for next time step 711 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) )712 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) )712 un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 713 vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 713 714 ! ! Save integrated transport for next computation 714 715 ub2_b(ji,jj) = zun_save … … 728 729 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 729 730 DO jk=1,jpkm1 730 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_ 2dt_b731 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_ 2dt_b731 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b 732 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b 732 733 END DO 733 734 ELSE … … 744 745 ! 745 746 DO jk=1,jpkm1 746 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_ 2dt_b747 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_ 2dt_b747 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 748 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 748 749 END DO 749 750 ! Save barotropic velocities not transport: … … 808 809 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 809 810 INTEGER, INTENT(inout) :: jpit ! cycle length 810 REAL(wp), DIMENSION(3*nn_ baro), INTENT(inout) :: zwgt1, & ! Primary weights811 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights 811 812 zwgt2 ! Secondary weights 812 813 … … 820 821 ! Set time index when averaged value is requested 821 822 IF (ll_fw) THEN 822 jic = nn_ baro823 jic = nn_e 823 824 ELSE 824 jic = 2 * nn_ baro825 jic = 2 * nn_e 825 826 ENDIF 826 827 … … 828 829 IF (ll_av) THEN 829 830 ! Define simple boxcar window for primary weights 830 ! (width = nn_ baro, centered around jic)831 ! (width = nn_e, centered around jic) 831 832 SELECT CASE ( nn_bt_flt ) 832 833 CASE( 0 ) ! No averaging … … 834 835 jpit = jic 835 836 836 CASE( 1 ) ! Boxcar, width = nn_ baro837 DO jn = 1, 3*nn_ baro838 za1 = ABS(float(jn-jic))/float(nn_ baro)837 CASE( 1 ) ! Boxcar, width = nn_e 838 DO jn = 1, 3*nn_e 839 za1 = ABS(float(jn-jic))/float(nn_e) 839 840 IF (za1 < 0.5_wp) THEN 840 841 zwgt1(jn) = 1._wp … … 843 844 ENDDO 844 845 845 CASE( 2 ) ! Boxcar, width = 2 * nn_ baro846 DO jn = 1, 3*nn_ baro847 za1 = ABS(float(jn-jic))/float(nn_ baro)846 CASE( 2 ) ! Boxcar, width = 2 * nn_e 847 DO jn = 1, 3*nn_e 848 za1 = ABS(float(jn-jic))/float(nn_e) 848 849 IF (za1 < 1._wp) THEN 849 850 zwgt1(jn) = 1._wp … … 889 890 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 890 891 ! ! --------------- 891 IF( ln_rstart .AND. ln_bt_fw .AND. ( neuler/=0) ) THEN !* Read the restart file892 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file 892 893 CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios ) 893 894 CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios ) … … 975 976 976 977 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 977 IF( ln_bt_auto ) nn_ baro = CEILING( rdt / rn_bt_cmax * zcmax)978 IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax) 978 979 979 r dtbt = rdt / REAL( nn_baro, wp )980 zcmax = zcmax * r dtbt980 rDt_e = rn_Dt / REAL( nn_e , wp ) 981 zcmax = zcmax * rDt_e 981 982 ! Print results 982 983 IF(lwp) WRITE(numout,*) … … 984 985 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 985 986 IF( ln_bt_auto ) THEN 986 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_ baro'987 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e ' 987 988 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 988 989 ELSE 989 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_ baro in namelist nn_baro = ', nn_baro990 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e 990 991 ENDIF 991 992 992 993 IF(ln_bt_av) THEN 993 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_ barotime steps is on '994 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on ' 994 995 ELSE 995 996 IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables ' … … 1011 1012 SELECT CASE ( nn_bt_flt ) 1012 1013 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1013 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_ baro'1014 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_ baro'1014 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e' 1015 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1015 1016 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' ) 1016 1017 END SELECT 1017 1018 ! 1018 1019 IF(lwp) WRITE(numout,*) ' ' 1019 IF(lwp) WRITE(numout,*) ' nn_ baro = ', nn_baro1020 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', r dtbt1020 IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e 1021 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e 1021 1022 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1022 1023 ! … … 1030 1031 ENDIF 1031 1032 IF( zcmax>0.9_wp ) THEN 1032 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_ baro!' )1033 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1033 1034 ENDIF 1034 1035 ! … … 1429 1430 ! 1430 1431 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1431 zztmp = -1._wp / r dtbt1432 zztmp = -1._wp / rDt_e 1432 1433 DO_2D_00_00 1433 1434 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & -
NEMO/trunk/src/OCE/DYN/dynzdf.F90
r12377 r12489 92 92 ENDIF 93 93 ENDIF 94 ! !* set time step95 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping)96 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog)97 ENDIF98 !99 94 ! !* explicit top/bottom drag case 100 95 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa)) … … 112 107 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 113 108 DO jk = 1, jpkm1 114 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r 2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk)115 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r 2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk)109 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + rDt * puu(:,:,jk,Krhs) ) * umask(:,:,jk) 110 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + rDt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk) 116 111 END DO 117 112 ELSE ! applied on thickness weighted velocity 118 113 DO jk = 1, jpkm1 119 114 puu(:,:,jk,Kaa) = ( e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) & 120 & + r 2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) ) / e3u(:,:,jk,Kaa) * umask(:,:,jk)115 & + rDt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) ) / e3u(:,:,jk,Kaa) * umask(:,:,jk) 121 116 pvv(:,:,jk,Kaa) = ( e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) & 122 & + r 2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk)117 & + rDt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk) 123 118 END DO 124 119 ENDIF … … 138 133 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 139 134 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 140 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r 2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua141 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r 2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va135 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 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 142 137 END_2D 143 138 IF( ln_isfcav ) THEN ! Ocean cavities (ISF) … … 147 142 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 148 143 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 149 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r 2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua150 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r 2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va144 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 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 151 146 END_2D 152 147 END IF … … 156 151 ! 157 152 ! !* Matrix construction 158 zdt = r 2dt * 0.5153 zdt = rDt * 0.5 159 154 IF( ln_zad_Aimp ) THEN !! 160 155 SELECT CASE( nldf_dyn ) … … 232 227 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 233 228 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 234 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r 2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua229 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 235 230 END_2D 236 231 IF ( ln_isfcav ) THEN ! top friction (always implicit) … … 239 234 iku = miku(ji,jj) ! ocean top level at u- and v-points 240 235 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 241 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r 2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua236 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 242 237 END_2D 243 238 END IF … … 265 260 DO_2D_00_00 266 261 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 267 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r 2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &268 & / ( ze3ua * r au0 ) * umask(ji,jj,1)262 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 263 & / ( ze3ua * rho0 ) * umask(ji,jj,1) 269 264 END_2D 270 265 DO_3D_00_00( 2, jpkm1 ) … … 282 277 ! 283 278 ! !* Matrix construction 284 zdt = r 2dt * 0.5279 zdt = rDt * 0.5 285 280 IF( ln_zad_Aimp ) THEN !! 286 281 SELECT CASE( nldf_dyn ) … … 357 352 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 358 353 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 359 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r 2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va354 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 360 355 END_2D 361 356 IF ( ln_isfcav ) THEN … … 363 358 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 364 359 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 365 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r 2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va360 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 366 361 END_2D 367 362 ENDIF … … 389 384 DO_2D_00_00 390 385 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 391 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r 2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &392 & / ( ze3va * r au0 ) * vmask(ji,jj,1)386 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 387 & / ( ze3va * rho0 ) * vmask(ji,jj,1) 393 388 END_2D 394 389 DO_3D_00_00( 2, jpkm1 ) … … 404 399 ! 405 400 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 406 ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r 2dt - ztrdu(:,:,:)407 ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r 2dt - ztrdv(:,:,:)401 ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:) 402 ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:) 408 403 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 409 404 DEALLOCATE( ztrdu, ztrdv ) -
NEMO/trunk/src/OCE/DYN/sshwzv.F90
r12377 r12489 75 75 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height 76 76 ! 77 INTEGER :: jk ! dummy loop indice78 REAL(wp) :: z 2dt, zcoef ! local scalars77 INTEGER :: jk ! dummy loop index 78 REAL(wp) :: zcoef ! local scalar 79 79 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 80 80 !!---------------------------------------------------------------------- … … 88 88 ENDIF 89 89 ! 90 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) 91 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 92 zcoef = 0.5_wp * r1_rau0 90 zcoef = 0.5_wp * r1_rho0 93 91 94 92 ! !------------------------------! … … 96 94 ! !------------------------------! 97 95 IF(ln_wd_il) THEN 98 CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv )96 CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv ) 99 97 ENDIF 100 98 … … 109 107 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 110 108 ! 111 pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)109 pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 112 110 ! 113 111 #if defined key_agrif … … 152 150 ! 153 151 INTEGER :: ji, jj, jk ! dummy loop indices 154 REAL(wp) :: z1_2dt ! local scalars155 152 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv 156 153 !!---------------------------------------------------------------------- … … 168 165 ! ! Now Vertical Velocity ! 169 166 ! !------------------------------! 170 z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog)171 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt172 167 ! 173 168 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases … … 187 182 ! computation of w 188 183 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk) & 189 & + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)184 & + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 190 185 END DO 191 186 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 … … 195 190 ! computation of w 196 191 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 197 & + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)192 & + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 198 193 END DO 199 194 ENDIF … … 227 222 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 228 223 !! from the filter, see Leclair and Madec 2010) and swap : 229 !! pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )230 !! - atfp * rdt * ( emp_b - emp ) / rau0224 !! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 225 !! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0 231 226 !! 232 227 !! ** action : - pssh(:,:,Kmm) time filtered … … 249 244 ENDIF 250 245 ! !== Euler time-stepping: no filter, just swap ==! 251 IF ( .NOT.( neuler == 0 .AND. kt == nit000) ) THEN ! Only do time filtering for leapfrog timesteps246 IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps 252 247 ! ! filtered "now" field 253 pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )248 pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 254 249 IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed 255 zcoef = atfp * rdt * r1_rau0250 zcoef = rn_atfp * rn_Dt * r1_rho0 256 251 pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) & 257 252 & - rnf_b(:,:) + rnf (:,:) & … … 260 255 261 256 ! ice sheet coupling 262 IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)257 IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:) 263 258 264 259 ENDIF … … 311 306 DO_3D_00_00( 1, jpkm1 ) 312 307 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 313 ! 2*r dt and not r2dt (for restartability)314 Cu_adv(ji,jj,jk) = 2._wp * r dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &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 ) ) & 315 310 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 316 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 ) ) & … … 324 319 DO_3D_00_00( 1, jpkm1 ) 325 320 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 326 ! 2*r dt and not r2dt (for restartability)327 Cu_adv(ji,jj,jk) = 2._wp * r dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &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 ) ) & 328 323 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 329 324 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & -
NEMO/trunk/src/OCE/DYN/wet_dry.F90
r12377 r12489 270 270 271 271 272 SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, r dtbt)272 SUBROUTINE wad_lmt_bt( zflxu, zflxv, sshn_e, zssh_frc, rDt_e ) 273 273 !!---------------------------------------------------------------------- 274 274 !! *** ROUTINE wad_lmt *** … … 280 280 !! ** Action : - calculate flux limiter and W/D flag 281 281 !!---------------------------------------------------------------------- 282 REAL(wp) , INTENT(in ) :: r dtbt! ocean time-step index282 REAL(wp) , INTENT(in ) :: rDt_e ! ocean time-step index 283 283 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zflxu, zflxv, sshn_e, zssh_frc 284 284 ! … … 299 299 zdepwd = 50._wp ! maximum depth that ocean cells can have W/D processes 300 300 ! 301 z2dt = r dtbt301 z2dt = rDt_e 302 302 ! 303 303 zflxp(:,:) = 0._wp
Note: See TracChangeset
for help on using the changeset viewer.
