Changeset 4178
- Timestamp:
- 2013-11-11T13:01:19+01:00 (10 years ago)
- Location:
- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO
- Files:
-
- 8 edited
-
LIM_SRC_2/limsbc_2.F90 (modified) (1 diff)
-
OPA_SRC/DYN/dynnxt.F90 (modified) (10 diffs)
-
OPA_SRC/DYN/dynspg.F90 (modified) (4 diffs)
-
OPA_SRC/DYN/dynspg_ts.F90 (modified) (19 diffs)
-
OPA_SRC/DYN/sshwzv.F90 (modified) (8 diffs)
-
OPA_SRC/SBC/updtide.F90 (modified) (3 diffs)
-
OPA_SRC/oce.F90 (modified) (2 diffs)
-
OPA_SRC/step.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
r3896 r4178 474 474 snwice_mass (:,:) = 0.e0 ! no mass exchanges 475 475 snwice_mass_b(:,:) = 0.e0 ! no mass exchanges 476 snwice_fmass (:,:) = 0.e0 ! no mass exchanges 476 477 ENDIF 477 478 IF( nn_ice_embd == 2 .AND. & ! full embedment (case 2) & no restart : -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r3865 r4178 17 17 !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module 18 18 !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 19 20 !!------------------------------------------------------------------------- 20 21 … … 42 43 USE wrk_nemo ! Memory Allocation 43 44 USE prtctl ! Print control 45 USE dynspg_ts ! Barotropic velocities 46 44 47 #if defined key_agrif 45 48 USE agrif_opa_interp … … 103 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 104 107 REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - - 108 REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva 105 109 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f 106 110 !!---------------------------------------------------------------------- … … 109 113 ! 110 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 115 IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva ) 111 116 ! 112 117 IF( kt == nit000 ) THEN … … 127 132 ! 128 133 #else 134 135 # if defined key_dynspg_exp 129 136 ! Next velocity : Leap-frog time stepping 130 137 ! ------------- … … 147 154 END DO 148 155 ENDIF 149 156 # endif 157 158 # if defined key_dynspg_ts 159 ! Ensure below that barotropic velocities match time splitting estimate 160 ! Compute actual transport and replace it with ts estimate at "after" time step 161 zua(:,:) = 0._wp 162 zva(:,:) = 0._wp 163 IF (lk_vvl) THEN 164 DO jk = 1, jpkm1 165 zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 167 END DO 168 DO jk = 1, jpkm1 169 ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_e(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_e(:,:) + va_b(:,:) ) * vmask(:,:,jk) 171 END DO 172 ELSE 173 DO jk = 1, jpkm1 174 zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 175 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 176 END DO 177 DO jk = 1, jpkm1 178 ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk) 179 va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk) 180 END DO 181 ENDIF 182 183 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 184 ! Remove advective velocity from "now velocities" 185 ! prior to asselin filtering 186 ! In the forward case, this is done below after asselin filtering 187 DO jk = 1, jpkm1 188 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 189 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 190 END DO 191 ENDIF 192 # endif 150 193 151 194 ! Update after velocity on domain lateral boundaries … … 194 237 vn(:,:,jk) = va(:,:,jk) 195 238 END DO 239 IF (lk_vvl) THEN 240 DO jk = 1, jpkm1 241 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 242 fse3u_b(:,:,jk) = fse3u_n(:,:,jk) 243 fse3v_b(:,:,jk) = fse3v_n(:,:,jk) 244 ENDDO 245 ENDIF 196 246 ELSE !* Leap-Frog : Asselin filter and swap 197 247 ! ! =============! … … 217 267 ! (used as a now filtered scale factor until the swap) 218 268 ! ---------------------------------------------------- 219 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 220 ! Add volume filter correction: compatibility with tracer advection scheme 221 ! => time filter + conservation correction (only at the first level) 222 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 269 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 270 ! Remove asselin filtering on thicknesses if forward time splitting 271 fse3t_b(:,:,:) = fse3t_n(:,:,:) 272 ELSE 273 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 274 ! Add volume filter correction: compatibility with tracer advection scheme 275 ! => time filter + conservation correction (only at the first level) 276 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 277 ! 278 ENDIF 223 279 ! 224 280 IF( ln_dynadv_vec ) THEN … … 276 332 ENDIF 277 333 ! 278 ENDIF 334 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 335 ! Remove asselin filtering of barotropic velocities if forward time splitting 336 ! note that we replace barotropic velocities by advective velocities 337 zua(:,:) = 0._wp 338 zva(:,:) = 0._wp 339 IF (lk_vvl) THEN 340 DO jk = 1, jpkm1 341 zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 342 zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 343 END DO 344 ELSE 345 DO jk = 1, jpkm1 346 zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 347 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 348 END DO 349 ENDIF 350 DO jk = 1, jpkm1 351 ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 352 vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 353 END DO 354 ENDIF 355 ! 356 ENDIF ! neuler =/0 279 357 280 358 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & … … 282 360 ! 283 361 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 362 IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva ) 284 363 ! 285 364 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt') -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r3957 r4178 22 22 USE dynspg_flt ! surface pressure gradient (dyn_spg_flt routine) 23 23 USE dynadv ! dynamics: vector invariant versus flux form 24 USE dynhpg, ONLY: ln_dynhpg_imp 24 25 USE sbctide 25 26 USE updtide … … 113 114 END DO 114 115 ! 115 IF( ln_apr_dyn ) THEN !== Atmospheric pressure gradient==!116 IF( ln_apr_dyn .AND. (.NOT. lk_dynspg_ts) ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==! 116 117 zg_2 = grav * 0.5 117 118 DO jj = 2, jpjm1 ! gradient of Patm using inverse barometer ssh … … 227 228 ENDIF 228 229 230 IF( lk_dynspg_ts ) CALL dyn_spg_ts_init( nit000 ) 231 ! (do it now, to set nn_baro, used to allocate some arrays later on) 229 232 ! ! allocate dyn_spg arrays 230 233 IF( lk_dynspg_ts ) THEN … … 266 269 ENDIF 267 270 268 ! ! Control of momentum formulation269 IF( lk_dynspg_ts .AND. l k_vvl) THEN270 IF( .NOT.ln_dynadv_vec ) CALL ctl_stop( 'Flux form not implemented for this free surface formulation' )271 ! ! Control of hydrostatic pressure choice 272 IF( lk_dynspg_ts .AND. ln_dynhpg_imp ) THEN 273 CALL ctl_stop( 'Semi-implicit hpg not compatible with time splitting' ) 271 274 ENDIF 272 275 ! -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3953 r4178 26 26 USE zdfbfr ! bottom friction 27 27 USE dynvor ! vorticity term 28 USE dynadv ! advection 28 29 USE obc_oce ! Lateral open boundary condition 29 30 USE obc_par ! open boundary condition parameters … … 34 35 USE bdydta ! open boundary condition data 35 36 USE bdydyn2d ! open boundary conditions on barotropic variables 37 USE bdytides ! 36 38 USE sbctide 37 39 USE updtide … … 49 51 PRIVATE 50 52 53 ! Potential namelist parameters below to be read in dyn_spg_ts_init 54 LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.TRUE. !: Forward integration of barotropic sub-stepping 55 LOGICAL, PRIVATE, PARAMETER :: ln_bt_av=.TRUE. !: Time averaging of barotropic variables 56 LOGICAL, PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE. !: Set number of iterations automatically 57 INTEGER, PRIVATE, PARAMETER :: nn_bt_flt=1 !: Filter choice 58 REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp !: Max. courant number (used if ln_bt_nn_auto=T) 59 ! End namelist parameters 60 61 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 62 REAL(wp),SAVE :: rdtbt ! Barotropic time step 63 64 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 65 wgtbtp1, & ! Primary weights used for time filtering of barotropic variables 66 wgtbtp2 ! Secondary weights used for time filtering of barotropic variables 67 ! 51 68 PUBLIC dyn_spg_ts ! routine called by step.F90 52 69 PUBLIC ts_rst ! routine called by istate.F90 53 70 PUBLIC dyn_spg_ts_alloc ! routine called by dynspg.F90 54 55 71 PUBLIC dyn_spg_ts_init ! " " " " 72 73 74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points 56 75 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 57 76 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 58 77 59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_b, vn_b ! now averaged velocity 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub_b, vb_b ! before averaged velocity 78 ! Would be convenient to have arrays below defined whatever the free surface option ? 79 ! These could be computed once for all at the beginning of the each baroclinic time step 80 ! and eventually swapped at the end 81 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ua_b, va_b ! after averaged velocities 82 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_b, vn_b ! now averaged velocities 83 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub_b, vb_b ! before averaged velocities 84 85 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv, vn_adv ! Advection vel. at "now" barocl. step 86 REAL(wp), PRIVATE,ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub2_b, vb2_b ! Advection vel. at "now-0.5" barocl. step 87 88 ! Arrays below are saved to allow testing of the "no time averaging" option 89 ! If this option is not retained, these could be replaced by temporary arrays 90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays 91 ubb_e, ub_e, & 92 vbb_e, vb_e 61 93 62 94 !! * Substitutions … … 74 106 !! *** routine dyn_spg_ts_alloc *** 75 107 !!---------------------------------------------------------------------- 76 ALLOCATE( ftnw (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) , & 77 & ftsw (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc ) 78 ! 108 INTEGER :: ierr(3) 109 !!---------------------------------------------------------------------- 110 ierr(:) = 0 111 112 ALLOCATE( ub_b(jpi,jpj) , vb_b(jpi,jpj) , & 113 & un_b(jpi,jpj) , vn_b(jpi,jpj) , & 114 & ua_b(jpi,jpj) , va_b(jpi,jpj) , & 115 & un_adv(jpi,jpj), vn_adv(jpi,jpj) , & 116 & sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 117 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 118 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , & 119 & wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), & 120 & zwz(jpi,jpj), STAT= ierr(1) ) 121 122 IF( ln_bt_fw ) ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), STAT=ierr(2) ) 123 124 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 125 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 126 127 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 128 79 129 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 80 130 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays') … … 113 163 INTEGER, INTENT(in) :: kt ! ocean time-step index 114 164 ! 165 LOGICAL :: ll_fw_start ! if true, forward integration 166 LOGICAL :: ll_init ! if true, special startup of 2d equations 115 167 INTEGER :: ji, jj, jk, jn ! dummy loop indices 116 INTEGER :: icycle ! local scalar 168 ! INTEGER :: icycle ! local scalar 169 INTEGER :: ioffset ! local scalar 117 170 INTEGER :: ikbu, ikbv ! local scalar 118 171 REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars 119 REAL(wp) :: z1_8, zx1, zy1 ! - - 120 REAL(wp) :: z1_4, zx2, zy2 ! - - 172 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 173 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 174 REAL(wp) :: za1, za2, za3 ! - - 121 175 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 122 176 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - … … 127 181 REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 128 182 REAL(wp), POINTER, DIMENSION(:,:) :: zhu_b, zhv_b 183 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zsshp2_e 129 184 !!---------------------------------------------------------------------- 130 185 ! … … 135 190 CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 136 191 CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b ) 137 ! 138 IF( kt == nit000 ) THEN !* initialisation 192 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e ) 193 ! 194 195 ! !* Local constant initialization 196 z1_12 = 1._wp / 12._wp 197 z1_8 = 0.125_wp ! coefficients for vorticity estimates 198 z1_4 = 0.25_wp 199 z1_2 = 0.5_wp 200 zraur = 1._wp / rau0 ! 1 / volumetric mass 201 ! 202 zhdiv(:,:) = 0._wp ! barotropic divergence 203 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 204 zv_sld = 0._wp ; zv_asp = 0._wp 205 206 IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction 207 z2dt_bf = rdt ! baroclinic time step (euler timestep) 208 ELSE 209 z2dt_bf = 2.0_wp * rdt ! baroclinic time step 210 ENDIF 211 z1_2dt_b = 1.0_wp / z2dt_bf ! reciprocal of baroclinic time step 212 ! 213 ll_init = ln_bt_av ! if no time averaging, then no specific restart 214 ll_fw_start = .FALSE. 215 ! 216 ! time offset in steps for bdy data update 217 IF (.NOT.ln_bt_fw) THEN ; ioffset=-2*nn_baro ; ELSE ; ioffset = 0 ; ENDIF 218 ! 219 IF( kt == nit000 ) THEN !* initialisation 139 220 ! 140 221 IF(lwp) WRITE(numout,*) … … 143 224 IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2*nn_baro 144 225 ! 145 CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b 226 IF (neuler==0) ll_init=.TRUE. 227 ! 228 IF (ln_bt_fw.OR.(neuler==0)) THEN 229 ll_fw_start=.TRUE. 230 ioffset = 0 231 ELSE 232 ll_fw_start=.FALSE. 233 ENDIF 234 ! 235 ! Set averaging weights and cycle length: 236 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 237 ! 238 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b 146 239 ! 147 240 ua_e (:,:) = un_b (:,:) … … 164 257 ! 165 258 ENDIF 166 167 ! !* Local constant initialization 168 z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! reciprocal of baroclinic time step 169 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! reciprocal of baroclinic 170 ! time step (euler timestep) 171 z1_8 = 0.125_wp ! coefficient for vorticity estimates 172 z1_4 = 0.25_wp 173 zraur = 1._wp / rau0 ! 1 / volumic mass 174 ! 175 zhdiv(:,:) = 0._wp ! barotropic divergence 176 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 177 zv_sld = 0._wp ; zv_asp = 0._wp 178 179 IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction 180 z2dt_bf = rdt 181 ELSE 182 z2dt_bf = 2.0_wp * rdt 259 ! 260 ! If forward start at previous time step, and centered integration, 261 ! then update averaging weights: 262 IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN 263 ll_fw_start=.FALSE. 264 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 183 265 ENDIF 184 266 … … 303 385 zcoef = -1._wp * z1_2dt_b 304 386 387 !--------> MERGED TO HERE 305 388 IF(ln_bfrimp) THEN 306 389 ! ! Remove the bottom stress trend from 3-D sea surface level gradient … … 358 441 ! ! ==================== ! 359 442 ! ! Initialisations ! 443 ! ! ==================== ! 444 ! Initialize barotropic variables: 445 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 446 sshn_e (:,:) = sshn (:,:) 447 zun_e (:,:) = un_b (:,:) 448 zvn_e (:,:) = vn_b (:,:) 449 ELSE ! CENTERED integration: start from BEFORE fields 450 sshn_e (:,:) = sshb (:,:) 451 zun_e (:,:) = ub_b (:,:) 452 zvn_e (:,:) = vb_b (:,:) 453 ENDIF 454 ! 455 ! Initialize depths: 456 IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 457 ! hu_e (:,:) = hu_0(:,:) + sshu_b(:,:) 458 ! hv_e (:,:) = hv_0(:,:) + sshv_b(:,:) 459 ! hur_e (:,:) = umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 460 ! hvr_e (:,:) = vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 461 hu_e (:,:) = hu (:,:) 462 hv_e (:,:) = hv (:,:) 463 hur_e (:,:) = hur (:,:) 464 hvr_e (:,:) = hvr (:,:) 465 ELSE 466 hu_e (:,:) = hu (:,:) 467 hv_e (:,:) = hv (:,:) 468 hur_e (:,:) = hur (:,:) 469 hvr_e (:,:) = hvr (:,:) 470 ENDIF 471 ! 472 IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5: 473 zhup2_e (:,:) = hu(:,:) 474 zhvp2_e (:,:) = hv(:,:) 475 ENDIF 476 ! 477 ! Initialize sums: 478 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) 479 va_b (:,:) = 0._wp 480 ssha (:,:) = 0._wp ! Sum for after averaged sea level 481 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop 482 zv_sum(:,:) = 0._wp 360 483 ! ! ==================== ! 361 icycle = 2 * nn_baro ! Number of barotropic sub time-step362 484 363 ! ! Start from NOW field364 hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points365 hv_e (:,:) = hv (:,:)366 hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points367 hvr_e (:,:) = hvr (:,:)368 !RBbug zsshb_e(:,:) = sshn (:,:)369 485 zsshb_e(:,:) = sshn_b(:,:) ! sea surface height (before and now) 370 sshn_e (:,:) = sshn (:,:)371 372 zun_e (:,:) = un_b (:,:) ! barotropic velocity (external)373 zvn_e (:,:) = vn_b (:,:)374 486 zub_e (:,:) = un_b (:,:) 375 487 zvb_e (:,:) = vn_b (:,:) … … 378 490 zv_sum (:,:) = vn_b (:,:) 379 491 zssh_sum(:,:) = sshn (:,:) 492 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) 493 va_b (:,:) = 0._wp 380 494 381 495 #if defined key_obc … … 396 510 ! !* Update the forcing (BDY and tides) 397 511 ! ! ------------------ 398 IF( lk_obc ) CALL obc_dta_bt( kt, jn ) 399 IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=1 ) 400 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, kbaro=nn_baro ) 401 402 ! !* after ssh_e 512 IF( lk_obc ) CALL obc_dta_bt( kt, jn ) 513 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(ioffset+1) ) 514 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=ioffset ) 515 ! 516 ! Set extrapolation coefficients for predictor step: 517 IF ((jn<3).AND.ll_init) THEN ! Forward 518 za1 = 1._wp 519 za2 = 0._wp 520 za3 = 0._wp 521 ELSE ! AB3-AM4 Coefficients: bet=0.281105 522 za1 = 1.781105_wp ! za1 = 3/2 + bet 523 za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet) 524 za3 = 0.281105_wp ! za3 = bet 525 ENDIF 526 527 ! Extrapolate barotropic velocities at step jit+0.5: 528 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 529 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 530 531 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 532 ! ! ------------------ 533 ! Extrapolate Sea Level at step jit+0.5: 534 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 535 ! 536 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 537 DO ji = 2, fs_jpim1 ! Vector opt. 538 zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 539 & * ( e1t(ji ,jj) * e2t(ji ,jj) * zsshp2_e(ji ,jj) & 540 & + e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 541 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 542 & * ( e1t(ji,jj ) * e2t(ji,jj ) * zsshp2_e(ji,jj ) & 543 & + e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 544 END DO 545 END DO 546 CALL lbc_lnk( zwx, 'U', 1._wp ) ; CALL lbc_lnk( zwy, 'V', 1._wp ) 547 ! 548 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 549 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 550 ENDIF 551 ! !* after ssh 403 552 ! ! ----------- 404 DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports 553 ! One should enforce volume conservation at open boundaries here 554 ! considering fluxes below: 555 ! 556 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 557 zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 558 DO jj = 2, jpjm1 405 559 DO ji = fs_2, fs_jpim1 ! vector opt. 406 560 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & … … 623 777 ENDIF 624 778 779 za1 = wgtbtp1(jn) 780 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities 781 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 782 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 783 ELSE ! Sum transports 784 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 785 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 786 ! ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 787 ! va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 788 ENDIF 789 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 790 625 791 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 626 792 ! ! ------------------ … … 664 830 zv_sum(:,:) = zcoef * zv_sum (:,:) 665 831 ! 832 ! Set advection velocity correction: 833 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 834 un_adv(:,:) = zu_sum(:,:) 835 vn_adv(:,:) = zv_sum(:,:) 836 ELSE 837 un_adv(:,:) = 0.5_wp * ( ub_b(:,:) + zu_sum(:,:)) 838 vn_adv(:,:) = 0.5_wp * ( vb_b(:,:) + zv_sum(:,:)) 839 END IF 666 840 ! !* update the general momentum trend 667 841 DO jk=1,jpkm1 … … 680 854 CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 681 855 CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b ) 856 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e ) 682 857 ! 683 858 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') … … 685 860 END SUBROUTINE dyn_spg_ts 686 861 862 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) 863 !!--------------------------------------------------------------------- 864 !! *** ROUTINE ts_wgt *** 865 !! 866 !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 867 !!---------------------------------------------------------------------- 868 LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true. 869 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 870 INTEGER, INTENT(inout) :: jpit ! cycle length 871 REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights 872 zwgt2 ! Secondary weights 873 874 INTEGER :: jic, jn, ji ! temporary integers 875 REAL(wp) :: za1, za2 876 !!---------------------------------------------------------------------- 877 878 zwgt1(:) = 0._wp 879 zwgt2(:) = 0._wp 880 881 ! Set time index when averaged value is requested 882 IF (ll_fw) THEN 883 jic = nn_baro 884 ELSE 885 jic = 2 * nn_baro 886 ENDIF 887 888 ! Set primary weights: 889 IF (ll_av) THEN 890 ! Define simple boxcar window for primary weights 891 ! (width = nn_baro, centered around jic) 892 SELECT CASE ( nn_bt_flt ) 893 CASE( 0 ) ! No averaging 894 zwgt1(jic) = 1._wp 895 jpit = jic 896 897 CASE( 1 ) ! Boxcar, width = nn_baro 898 DO jn = 1, 3*nn_baro 899 za1 = ABS(float(jn-jic))/float(nn_baro) 900 IF (za1 < 0.5_wp) THEN 901 zwgt1(jn) = 1._wp 902 jpit = jn 903 ENDIF 904 ENDDO 905 906 CASE( 2 ) ! Boxcar, width = 2 * nn_baro 907 DO jn = 1, 3*nn_baro 908 za1 = ABS(float(jn-jic))/float(nn_baro) 909 IF (za1 < 1._wp) THEN 910 zwgt1(jn) = 1._wp 911 jpit = jn 912 ENDIF 913 ENDDO 914 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 915 END SELECT 916 917 ELSE ! No time averaging 918 zwgt1(jic) = 1._wp 919 jpit = jic 920 ENDIF 921 922 ! Set secondary weights 923 DO jn = 1, jpit 924 DO ji = jn, jpit 925 zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 926 END DO 927 END DO 928 929 ! Normalize weigths: 930 za1 = 1._wp / SUM(zwgt1(1:jpit)) 931 za2 = 1._wp / SUM(zwgt2(1:jpit)) 932 DO jn = 1, jpit 933 zwgt1(jn) = zwgt1(jn) * za1 934 zwgt2(jn) = zwgt2(jn) * za2 935 END DO 936 ! 937 END SUBROUTINE ts_wgt 687 938 688 939 SUBROUTINE ts_rst( kt, cdrw ) … … 809 1060 END SUBROUTINE ts_rst 810 1061 1062 SUBROUTINE dyn_spg_ts_init( kt ) 1063 !!--------------------------------------------------------------------- 1064 !! *** ROUTINE dyn_spg_ts_init *** 1065 !! 1066 !! ** Purpose : Set time splitting options 1067 !!---------------------------------------------------------------------- 1068 INTEGER , INTENT(in) :: kt ! ocean time-step 1069 ! 1070 INTEGER :: ji ,jj 1071 REAL(wp) :: zxr2, zyr2, zcmax 1072 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1073 !! 1074 ! NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, & 1075 ! & nn_baro, rn_bt_cmax, nn_bt_flt 1076 !!---------------------------------------------------------------------- 1077 ! REWIND( numnam ) !* Namelist namsplit: split-explicit free surface 1078 ! READ ( numnam, namsplit ) 1079 ! ! Max courant number for ext. grav. waves 1080 ! 1081 CALL wrk_alloc( jpi, jpj, zcu ) 1082 ! 1083 IF (lk_vvl) THEN 1084 DO jj = 1, jpj 1085 DO ji =1, jpi 1086 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1087 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1088 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1089 END DO 1090 END DO 1091 ELSE 1092 DO jj = 1, jpj 1093 DO ji =1, jpi 1094 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1095 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1096 zcu(ji,jj) = sqrt( grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1097 END DO 1098 END DO 1099 ENDIF 1100 1101 zcmax = MAXVAL(zcu(:,:)) 1102 IF( lk_mpp ) CALL mpp_max( zcmax ) 1103 1104 ! Estimate number of iterations to satisfy a max courant number=0.8 1105 IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1106 1107 rdtbt = rdt / FLOAT(nn_baro) 1108 zcmax = zcmax * rdtbt 1109 ! Print results 1110 IF(lwp) WRITE(numout,*) 1111 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1112 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1113 IF( ln_bt_nn_auto ) THEN 1114 IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro ' 1115 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1116 ELSE 1117 IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist ' 1118 ENDIF 1119 IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro 1120 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt 1121 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1122 1123 IF(ln_bt_av) THEN 1124 IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on ' 1125 ELSE 1126 IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables ' 1127 ENDIF 1128 ! 1129 ! 1130 IF(ln_bt_fw) THEN 1131 IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables ' 1132 ELSE 1133 IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables ' 1134 ENDIF 1135 ! 1136 IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 1137 SELECT CASE ( nn_bt_flt ) 1138 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1139 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1140 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1141 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1142 END SELECT 1143 ! 1144 IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN 1145 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 1146 ENDIF 1147 IF ( zcmax>0.9_wp ) THEN 1148 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1149 ENDIF 1150 ! 1151 CALL wrk_dealloc( jpi, jpj, zcu ) 1152 ! 1153 END SUBROUTINE dyn_spg_ts_init 1154 811 1155 #else 812 1156 !!---------------------------------------------------------------------- 813 1157 !! Default case : Empty module No standart free surface cst volume 814 1158 !!---------------------------------------------------------------------- 1159 1160 USE par_kind 1161 LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping 1162 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_b, vn_b ! declaration needed for compilation when "key_dynspg_ts" is not defined. Arrays never used or allocated 1163 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv, vn_adv ! declaration needed for compilation when "key_dynspg_ts" is not defined. Arrays never used or allocated 1164 815 1165 CONTAINS 816 1166 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function … … 826 1176 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw 827 1177 END SUBROUTINE ts_rst 1178 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine 1179 INTEGER , INTENT(in) :: kt ! ocean time-step 1180 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt 1181 END SUBROUTINE dyn_spg_ts_init 828 1182 #endif 829 1183 -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r3896 r4178 33 33 USE diaar5, ONLY: lk_diaar5 34 34 USE iom 35 USE sbcrnf, ONLY: h_rnf, nk_rnf ! River runoff 35 USE sbcrnf, ONLY: h_rnf, nk_rnf, sbc_rnf_div ! River runoff 36 USE dynspg_ts, ONLY: un_b, vn_b, un_adv, vn_adv 37 USE dynspg_oce, ONLY: lk_dynspg_ts 36 38 #if defined key_agrif 37 39 USE agrif_opa_update … … 48 50 49 51 PUBLIC ssh_nxt ! called by step.F90 50 PUBLIC wzv ! called by step.F90 52 PUBLIC wzv_1 ! called by step.F90 53 PUBLIC wzv_2 ! called by step.F90 51 54 PUBLIC ssh_swp ! called by step.F90 52 55 … … 149 152 150 153 151 SUBROUTINE wzv ( kt )152 !!---------------------------------------------------------------------- 153 !! *** ROUTINE wzv ***154 SUBROUTINE wzv_1( kt ) 155 !!---------------------------------------------------------------------- 156 !! *** ROUTINE wzv_1 *** 154 157 !! 155 158 !! ** Purpose : compute the now vertical velocity … … 166 169 ! 167 170 INTEGER, INTENT(in) :: kt ! time step 168 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 169 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 171 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhdiv 170 172 ! 171 173 INTEGER :: ji, jj, jk ! dummy loop indices … … 173 175 !!---------------------------------------------------------------------- 174 176 175 IF( nn_timing == 1 ) CALL timing_start('wzv') 176 ! 177 CALL wrk_alloc( jpi, jpj, z2d ) 178 CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 177 IF( nn_timing == 1 ) CALL timing_start('wzv_1') 178 ! 179 CALL wrk_alloc( jpi, jpj, jpk, zhdiv ) 179 180 ! 180 181 IF( kt == nit000 ) THEN 181 182 ! 182 183 IF(lwp) WRITE(numout,*) 183 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '184 IF(lwp) WRITE(numout,*) '~~~ '184 IF(lwp) WRITE(numout,*) 'wzv_1 : now vertical velocity ' 185 IF(lwp) WRITE(numout,*) '~~~~~ ' 185 186 ! 186 187 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) … … 224 225 #endif 225 226 227 ! 228 CALL wrk_dealloc( jpi, jpj, jpk, zhdiv ) 229 ! 230 IF( nn_timing == 1 ) CALL timing_stop('wzv_1') 231 232 233 END SUBROUTINE wzv_1 234 235 SUBROUTINE wzv_2( kt ) 236 !!---------------------------------------------------------------------- 237 !! *** ROUTINE wzv_2 *** 238 !! 239 !! ** Purpose : compute the now vertical velocity 240 !! 241 !! ** Method : - Using the incompressibility hypothesis, the vertical 242 !! velocity is computed by integrating the horizontal divergence 243 !! from the bottom to the surface minus the scale factor evolution. 244 !! The boundary conditions are w=0 at the bottom (no flux) and. 245 !! 246 !! ** action : wn : now vertical velocity 247 !! 248 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 249 !!---------------------------------------------------------------------- 250 ! 251 INTEGER, INTENT(in) :: kt ! time step 252 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 253 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 254 ! 255 INTEGER :: ji, jj, jk ! dummy loop indices 256 REAL(wp) :: z1_2dt ! local scalars 257 !!---------------------------------------------------------------------- 258 259 IF( nn_timing == 1 ) CALL timing_start('wzv_2') 260 ! 261 CALL wrk_alloc( jpi, jpj, jpk, z3d, zhdiv ) 262 ! 263 IF( kt == nit000 ) THEN 264 ! 265 IF(lwp) WRITE(numout,*) 266 IF(lwp) WRITE(numout,*) 'wzv_2 : now vertical velocity ' 267 IF(lwp) WRITE(numout,*) '~~~~~ ' 268 ! 269 ENDIF 270 ! !------------------------------! 271 ! ! Now Vertical Velocity ! 272 ! !------------------------------! 273 z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog) 274 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt 275 276 IF( lk_dynspg_ts ) THEN 277 ! At this stage: 278 ! 1) vertical scale factors have been updated. 279 ! 2) Time averaged barotropic velocities around step kt+1 (ua_b, va_b) as well as 280 ! "advective" barotropic velocities (un_adv, vn_adv) at step kt in agreement 281 ! with continuity equation are available. 282 ! 3) 3d velocities (ua, va) hold ua_b, va_b issued from time splitting as barotropic components. 283 ! Below => Update "now" velocities, divergence, then vertical velocity 284 ! NB: hdivn is NOT updated such that hdivb is not updated also. This means that horizontal 285 ! momentum diffusion is still performed on Instantaneous barotropic arrays. I experencied 286 ! some issues with UBS with the current method. See "divup" comments in development branch to 287 ! update the divergence fully if necessary (2013/dev_r3867_MERCATOR1_DYN). 288 ! 289 DO jk = 1, jpkm1 290 ! Correct velocities: 291 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 292 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 293 ! 294 ! Compute divergence: 295 DO jj = 2, jpjm1 296 DO ji = fs_2, fs_jpim1 ! vector opt. 297 z3d(ji,jj,jk) = & 298 ( e2u(ji,jj)*fse3u(ji,jj,jk) * un(ji,jj,jk) - e2u(ji-1,jj)*fse3u(ji-1,jj,jk) * un(ji-1,jj,jk) & 299 + e1v(ji,jj)*fse3v(ji,jj,jk) * vn(ji,jj,jk) - e1v(ji,jj-1)*fse3v(ji,jj-1,jk) * vn(ji,jj-1,jk) ) & 300 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 301 END DO 302 END DO 303 END DO 304 ! 305 IF( ln_rnf ) CALL sbc_rnf_div( z3d ) ! runoffs (update divergence) 306 ELSE 307 z3d(:,:,:) = hdivn(:,:,:) 308 ENDIF 309 ! 310 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 311 DO jk = 1, jpkm1 312 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 313 ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 314 DO jj = 2, jpjm1 315 DO ji = fs_2, fs_jpim1 ! vector opt. 316 zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 317 END DO 318 END DO 319 END DO 320 CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 321 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 322 ! ! Same question holds for hdivn. Perhaps just for security 323 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 324 ! computation of w 325 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * z3d(:,:,jk) + zhdiv(:,:,jk) & 326 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 327 END DO 328 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 329 ELSE ! z_star and linear free surface cases 330 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 331 ! computation of w 332 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * z3d(:,:,jk) & 333 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 334 END DO 335 ENDIF 336 337 #if defined key_bdy 338 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 339 #endif 340 226 341 ! !------------------------------! 227 342 ! ! outputs ! … … 229 344 CALL iom_put( "woce", wn ) ! vertical velocity 230 345 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 346 CALL wrk_alloc( jpi, jpj, z2d ) 231 347 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 232 348 z2d(:,:) = rau0 * e12t(:,:) … … 236 352 CALL iom_put( "w_masstr" , z3d ) 237 353 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 238 ENDIF239 !240 CALL wrk_dealloc( jpi, jpj, z2d )354 CALL wrk_dealloc( jpi, jpj, z2d ) 355 ENDIF 356 ! 241 357 CALL wrk_dealloc( jpi, jpj, jpk, z3d, zhdiv ) 242 358 ! 243 IF( nn_timing == 1 ) CALL timing_stop('wzv') 244 245 246 END SUBROUTINE wzv 247 359 IF( nn_timing == 1 ) CALL timing_stop('wzv_2') 360 361 362 END SUBROUTINE wzv_2 248 363 249 364 SUBROUTINE ssh_swp( kt ) -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/SBC/updtide.F90
r3953 r4178 31 31 CONTAINS 32 32 33 SUBROUTINE upd_tide( kt, kit, kbaro )33 SUBROUTINE upd_tide( kt, kit, kbaro, koffset ) 34 34 !!---------------------------------------------------------------------- 35 35 !! *** ROUTINE upd_tide *** … … 44 44 INTEGER, INTENT(in), OPTIONAL :: kit ! external mode sub-time-step index (lk_dynspg_ts=T only) 45 45 INTEGER, INTENT(in), OPTIONAL :: kbaro ! number of sub-time-step (lk_dynspg_ts=T only) 46 INTEGER, INTENT(in), OPTIONAL :: koffset ! time offset in number 47 ! of sub-time-steps (lk_dynspg_ts=T only) 46 48 ! 49 INTEGER :: joffset ! local integer 47 50 INTEGER :: ji, jj, jk ! dummy loop indices 48 51 REAL(wp) :: zt, zramp ! local scalar … … 52 55 ! ! tide pulsation at model time step (or sub-time-step) 53 56 zt = ( kt - kt_tide ) * rdt 54 IF( PRESENT( kit ) .AND. PRESENT( kbaro ) ) zt = zt + kit * rdt / REAL( kbaro, wp ) 57 ! 58 joffset = 0 59 IF( PRESENT( koffset ) ) joffset = koffset 60 ! 61 IF( PRESENT( kit ) .AND. PRESENT( kbaro ) ) THEN 62 zt = zt + ( kit + 0.5_wp * ( joffset - 1 ) ) * rdt / REAL( kbaro, wp ) 63 ELSE 64 zt = zt + joffset * rdt 65 ENDIF 66 ! 55 67 zwt(:) = omega_tide(:) * zt 56 68 -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/oce.F90
r3865 r4178 22 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ub , un , ua !: i-horizontal velocity [m/s] 23 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vb , vn , va !: j-horizontal velocity [m/s] 24 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ua_sv, va_sv !: Saved trends (time spliting) [m/s2] 24 25 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wn !: vertical velocity [m/s] 25 26 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rotb , rotn !: relative vorticity [s-1] … … 68 69 ALLOCATE( ub (jpi,jpj,jpk) , un (jpi,jpj,jpk) , ua(jpi,jpj,jpk) , & 69 70 & vb (jpi,jpj,jpk) , vn (jpi,jpj,jpk) , va(jpi,jpj,jpk) , & 71 & ua_sv(jpi,jpj,jpk) , va_sv(jpi,jpj,jpk) , & 70 72 & wn (jpi,jpj,jpk) , & 71 73 & rotb (jpi,jpj,jpk) , rotn (jpi,jpj,jpk) , & -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/step.F90
r3953 r4178 99 99 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 100 100 CALL ssh_nxt ( kstp ) ! after ssh 101 IF( lk_dynspg_ts ) THEN 102 CALL wzv_1 ( kstp ) ! now cross-level velocity 103 ! In case the time splitting case, update almost all momentum trends here: 104 ! Note that the computation of vertical velocity above, hence "after" sea level 105 ! is necessary to compute momentum advection for the rhs of barotropic loop: 106 CALL eos ( tsn, rhd, rhop ) ! now in situ density for hpg computation 107 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 108 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 109 CALL zdf_bfr( kstp ) ! bottom friction (if quadratic) 110 111 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 112 va(:,:,:) = 0.e0 113 IF( ln_asmiau .AND. & 114 & ln_dyninc ) CALL dyn_asm_inc ( kstp ) ! apply dynamics assimilation increment 115 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! subtract Neptune velocities (simplified) 116 IF( lk_bdy ) CALL bdy_dyn3d_dmp( kstp ) ! bdy damping trends 117 CALL dyn_adv ( kstp ) ! advection (vector or flux form) 118 CALL dyn_vor ( kstp ) ! vorticity term including Coriolis 119 CALL dyn_ldf ( kstp ) ! lateral mixing 120 IF( ln_neptsimp ) CALL dyn_nept_cor ( kstp ) ! add Neptune velocities (simplified) 121 #if defined key_agrif 122 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momentum sponge 123 #endif 124 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure 125 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 126 127 ua_sv(:,:,:) = ua(:,:,:) ! save next velocities (not trends !) 128 va_sv(:,:,:) = va(:,:,:) 129 ENDIF 101 130 IF( lk_vvl ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 102 CALL wzv ( kstp ) ! now cross-level velocity131 CALL wzv_2 ( kstp ) ! now cross-level velocity (original) 103 132 104 133 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 109 138 ! 110 139 ! VERTICAL PHYSICS 111 140 IF( .NOT. lk_dynspg_ts ) CALL zdf_bfr( kstp ) ! bottom friction 112 141 113 142 ! ! Vertical eddy viscosity and diffusivity coefficients … … 208 237 209 238 ELSE ! centered hpg (eos then time stepping) 210 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 211 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 239 IF ( .NOT. lk_dynspg_ts ) THEN ! eos already called in time-split case 240 CALL eos ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation 241 IF( ln_zps ) CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv, & ! zps: now hor. derivative 212 242 & rhd, gru , grv ) ! of t, s, rd at the last ocean level 243 ENDIF 213 244 IF( ln_zdfnpc ) CALL tra_npc( kstp ) ! update after fields by non-penetrative convection 214 245 CALL tra_nxt( kstp ) ! tracer fields at next time step … … 218 249 ! Dynamics (tsa used as workspace) 219 250 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 251 IF( lk_dynspg_ts ) THEN 252 ! revert to previously computed tendencies: 253 ! (not using ua, va as temporary arrays during tracers' update could avoid that) 254 ua(:,:,:) = ua_sv(:,:,:) 255 va(:,:,:) = va_sv(:,:,:) 256 CALL dyn_bfr( kstp ) ! bottom friction 257 CALL dyn_zdf( kstp ) ! vertical diffusion 258 ELSE 220 259 ua(:,:,:) = 0.e0 ! set dynamics trends to zero 221 260 va(:,:,:) = 0.e0 222 261 223 IF( ln_asmiau .AND. &224 & ln_dyninc )CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment225 IF( ln_bkgwri )CALL asm_bkg_wri( kstp ) ! output background fields226 IF( ln_neptsimp )CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified)227 IF( lk_bdy )CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends262 IF( ln_asmiau .AND. & 263 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 264 IF( ln_bkgwri ) CALL asm_bkg_wri( kstp ) ! output background fields 265 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 266 IF( lk_bdy ) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 228 267 CALL dyn_adv( kstp ) ! advection (vector or flux form) 229 268 CALL dyn_vor( kstp ) ! vorticity term including Coriolis 230 269 CALL dyn_ldf( kstp ) ! lateral mixing 231 IF( ln_neptsimp )CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified)232 #if defined key_agrif 233 IF(.NOT. Agrif_Root())CALL Agrif_Sponge_dyn ! momemtum sponge270 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! add Neptune velocities (simplified) 271 #if defined key_agrif 272 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn ! momemtum sponge 234 273 #endif 235 274 CALL dyn_hpg( kstp ) ! horizontal gradient of Hydrostatic pressure … … 237 276 CALL dyn_zdf( kstp ) ! vertical diffusion 238 277 CALL dyn_spg( kstp, indic ) ! surface pressure gradient 278 ENDIF 239 279 CALL dyn_nxt( kstp ) ! lateral velocity at next time step 240 280
Note: See TracChangeset
for help on using the changeset viewer.
