- Timestamp:
- 2011-10-12T14:28:01+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2739_LOCEAN8_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r2715 r2905 8 8 !! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 10 11 !!---------------------------------------------------------------------- 11 12 … … 43 44 PRIVATE 44 45 45 PUBLIC ssh_wzv ! called by step.F9046 46 PUBLIC ssh_nxt ! called by step.F90 47 PUBLIC wzv ! called by step.F90 48 PUBLIC ssh_swp ! called by step.F90 47 49 48 50 !! * Substitutions … … 56 58 CONTAINS 57 59 58 SUBROUTINE ssh_ wzv( kt )59 !!---------------------------------------------------------------------- 60 !! *** ROUTINE ssh_ wzv***60 SUBROUTINE ssh_nxt( kt ) 61 !!---------------------------------------------------------------------- 62 !! *** ROUTINE ssh_nxt *** 61 63 !! 62 !! ** Purpose : compute the after ssh (ssha), the now vertical velocity 63 !! and update the now vertical coordinate (lk_vvl=T). 64 !! 65 !! ** Method : - Using the incompressibility hypothesis, the vertical 66 !! velocity is computed by integrating the horizontal divergence 67 !! from the bottom to the surface minus the scale factor evolution. 68 !! The boundary conditions are w=0 at the bottom (no flux) and. 64 !! ** Purpose : compute the after ssh (ssha) 65 !! 66 !! ** Method : - Using the incompressibility hypothesis, the ssh increment 67 !! is computed by integrating the horizontal divergence and multiply by 68 !! by the time step. 69 69 !! 70 70 !! ** action : ssha : after sea surface height 71 !! wn : now vertical velocity72 !! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T)73 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points74 71 !! 75 72 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 76 73 !!---------------------------------------------------------------------- 77 74 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 78 USE oce , ONLY: z3d => ta ! ta used as 3D workspace 79 USE wrk_nemo, ONLY: zhdiv => wrk_2d_1 , z2d => wrk_2d_2 ! 2D workspace 80 ! 81 INTEGER, INTENT(in) :: kt ! time step 82 ! 83 INTEGER :: ji, jj, jk ! dummy loop indices 84 REAL(wp) :: zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0 ! local scalars 85 !!---------------------------------------------------------------------- 86 87 IF( wrk_in_use(2, 1,2) ) THEN 88 CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable') ; RETURN 75 USE wrk_nemo, ONLY: zhdiv => wrk_2d_1 ! 2D workspace 76 ! 77 INTEGER, INTENT(in) :: kt ! time step 78 ! 79 INTEGER :: jk ! dummy loop indice 80 REAL(wp) :: z2dt, z1_rau0 ! local scalars 81 !!---------------------------------------------------------------------- 82 83 IF( wrk_in_use(2, 1) ) THEN 84 CALL ctl_stop('ssh_nxt: requested workspace arrays unavailable') ; RETURN 89 85 ENDIF 90 86 … … 92 88 ! 93 89 IF(lwp) WRITE(numout,*) 94 IF(lwp) WRITE(numout,*) 'ssh_ wzv : after sea surface height and now vertical velocity'90 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 95 91 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 96 92 ! 97 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)98 !99 IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only)100 DO jj = 1, jpjm1101 DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible102 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )103 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )104 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)105 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) &106 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )107 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) &108 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )109 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) &110 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )111 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) &112 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )113 END DO114 END DO115 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. )116 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. )117 DO jj = 1, jpjm1118 DO ji = 1, jpim1 ! NO Vector Opt.119 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) &120 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) &121 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) &122 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )123 END DO124 END DO125 CALL lbc_lnk( sshf_n, 'F', 1. )126 ENDIF127 !128 ENDIF129 130 ! !------------------------------------------!131 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case)132 ! !------------------------------------------!133 DO jk = 1, jpkm1134 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays135 fsdepw(:,:,jk) = fsdepw_n(:,:,jk)136 fsde3w(:,:,jk) = fsde3w_n(:,:,jk)137 !138 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays139 fse3u (:,:,jk) = fse3u_n (:,:,jk)140 fse3v (:,:,jk) = fse3v_n (:,:,jk)141 fse3f (:,:,jk) = fse3f_n (:,:,jk)142 fse3w (:,:,jk) = fse3w_n (:,:,jk)143 fse3uw(:,:,jk) = fse3uw_n(:,:,jk)144 fse3vw(:,:,jk) = fse3vw_n(:,:,jk)145 END DO146 !147 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points)148 hv(:,:) = hv_0(:,:) + sshv_n(:,:)149 ! ! now masked inverse of the ocean depth (at u- and v-points)150 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )151 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )152 !153 93 ENDIF 154 94 ! … … 184 124 CALL lbc_lnk( ssha, 'T', 1. ) 185 125 #endif 186 187 ! ! Sea Surface Height at u-,v- and f-points (vvl case only)188 IF( lk_vvl ) THEN ! (required only in key_vvl case)189 DO jj = 1, jpjm1190 DO ji = 1, jpim1 ! NO Vector Opt.191 sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) &192 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) &193 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )194 sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) &195 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) &196 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )197 END DO198 END DO199 CALL lbc_lnk( sshu_a, 'U', 1. ) ; CALL lbc_lnk( sshv_a, 'V', 1. ) ! Boundaries conditions200 ENDIF201 202 126 #if defined key_asminc 203 127 ! ! Include the IAU weighted SSH increment … … 209 133 210 134 ! !------------------------------! 135 ! ! outputs ! 136 ! !------------------------------! 137 CALL iom_put( "ssh" , sshn ) ! sea surface height 138 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 139 ! 140 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 141 ! 142 IF( wrk_not_released(2, 1) ) CALL ctl_stop('ssh_nxt: failed to release workspace arrays') 143 ! 144 END SUBROUTINE ssh_nxt 145 146 147 SUBROUTINE wzv( kt ) 148 !!---------------------------------------------------------------------- 149 !! *** ROUTINE wzv *** 150 !! 151 !! ** Purpose : compute the now vertical velocity 152 !! 153 !! ** Method : - Using the incompressibility hypothesis, the vertical 154 !! velocity is computed by integrating the horizontal divergence 155 !! from the bottom to the surface minus the scale factor evolution. 156 !! The boundary conditions are w=0 at the bottom (no flux) and. 157 !! 158 !! ** action : wn : now vertical velocity 159 !! 160 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 161 !!---------------------------------------------------------------------- 162 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 163 USE oce , ONLY: z3d => ta ! ta used as 3D workspace 164 USE wrk_nemo, ONLY: zhdiv => wrk_3d_1 , z2d => wrk_2d_1 ! 2D workspace 165 ! 166 INTEGER, INTENT(in) :: kt ! time step 167 ! 168 INTEGER :: ji, jj, jk ! dummy loop indices 169 REAL(wp) :: z1_2dt ! local scalars 170 !!---------------------------------------------------------------------- 171 172 IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 1) ) THEN 173 CALL ctl_stop('wzv: requested workspace arrays unavailable') ; RETURN 174 ENDIF 175 176 IF( kt == nit000 ) THEN 177 ! 178 IF(lwp) WRITE(numout,*) 179 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' 180 IF(lwp) WRITE(numout,*) '~~~ ' 181 ! 182 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 183 ! 184 ENDIF 185 ! !------------------------------! 211 186 ! ! Now Vertical Velocity ! 212 187 ! !------------------------------! 213 z1_2dt = 1.e0 / z2dt 214 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 215 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 216 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 217 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 218 & * tmask(:,:,jk) * z1_2dt 188 z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog) 189 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt 190 ! 191 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases 192 DO jk = 1, jpkm1 193 ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t) 194 ! - ML - note: computation allready done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way) 195 DO jj = 2, jpjm1 196 DO ji = fs_2, fs_jpim1 ! vector opt. 197 zhdiv(ji,jj,jk) = e12t_1(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 198 END DO 199 END DO 200 END DO 201 CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 202 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 203 ! ! Same question holds for hdivn. Perhaps just for security 204 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 205 ! computation of w 206 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) & 207 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 208 END DO 209 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 210 ELSE ! z_star and linear free surface cases 211 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 212 ! computation of w 213 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) & 214 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk) 215 END DO 216 ENDIF 217 219 218 #if defined key_bdy 220 219 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 221 220 #endif 222 END DO223 221 224 222 ! !------------------------------! 225 223 ! ! outputs ! 226 224 ! !------------------------------! 227 CALL iom_put( "woce", wn ) ! vertical velocity 228 CALL iom_put( "ssh" , sshn ) ! sea surface height 229 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 225 CALL iom_put( "woce", wn ) ! vertical velocity 230 226 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 231 227 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. … … 237 233 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 238 234 ENDIF 239 ! 240 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 241 ! 242 IF( wrk_not_released(2, 1,2) ) CALL ctl_stop('ssh_wzv: failed to release workspace arrays') 243 ! 244 END SUBROUTINE ssh_wzv 245 246 247 SUBROUTINE ssh_nxt( kt ) 235 236 IF( wrk_not_released(2, 1) .OR. wrk_not_released(3, 1) ) CALL ctl_stop('wzv: failed to release workspace arrays') 237 238 END SUBROUTINE wzv 239 240 241 SUBROUTINE ssh_swp( kt ) 248 242 !!---------------------------------------------------------------------- 249 243 !! *** ROUTINE ssh_nxt *** … … 251 245 !! ** Purpose : achieve the sea surface height time stepping by 252 246 !! applying Asselin time filter and swapping the arrays 253 !! ssha already computed in ssh_ wzv247 !! ssha already computed in ssh_nxt 254 248 !! 255 249 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing … … 266 260 INTEGER, INTENT(in) :: kt ! ocean time-step index 267 261 !! 268 INTEGER :: ji, jj ! dummy loop indices 269 REAL(wp) :: zec ! temporary scalar 262 REAL(wp) :: zec ! temporary scalar 270 263 !!---------------------------------------------------------------------- 271 264 272 265 IF( kt == nit000 ) THEN 273 266 IF(lwp) WRITE(numout,*) 274 IF(lwp) WRITE(numout,*) 'ssh_ nxt : next sea surface height (Asselin time filter + swap)'267 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height' 275 268 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 276 269 ENDIF 277 270 278 ! !--------------------------! 279 IF( lk_vvl ) THEN ! Variable volume levels ! (ssh at t-, u-, v, f-points) 280 ! !--------------------------! 281 ! 282 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 283 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 284 sshu_n(:,:) = sshu_a(:,:) 285 sshv_n(:,:) = sshv_a(:,:) 286 DO jj = 1, jpjm1 ! ssh now at f-point 287 DO ji = 1, jpim1 ! NO Vector Opt. 288 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 289 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 290 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 291 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 292 END DO 293 END DO 294 CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions 295 ! 296 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 271 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 272 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 273 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 274 IF( lk_vvl ) THEN ! before <-- now filtered 297 275 zec = atfp * rdt / rau0 298 DO jj = 1, jpj 299 DO ji = 1, jpi ! before <-- now filtered 300 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 301 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 302 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 303 sshu_n(ji,jj) = sshu_a(ji,jj) 304 sshv_n(ji,jj) = sshv_a(ji,jj) 305 END DO 306 END DO 307 DO jj = 1, jpjm1 ! ssh now at f-point 308 DO ji = 1, jpim1 ! NO Vector Opt. 309 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 310 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 311 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 312 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 313 END DO 314 END DO 315 CALL lbc_lnk( sshf_n, 'F', 1. ) ! Boundaries conditions 316 ! 317 DO jj = 1, jpjm1 ! ssh before at u- & v-points 318 DO ji = 1, jpim1 ! NO Vector Opt. 319 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 320 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 321 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 322 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 323 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 324 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 325 END DO 326 END DO 327 CALL lbc_lnk( sshu_b, 'U', 1. ) 328 CALL lbc_lnk( sshv_b, 'V', 1. ) ! Boundaries conditions 329 ! 276 sshb (:,:) = sshn (:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) & 277 & - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 278 ELSE 279 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 330 280 ENDIF 331 ! !--------------------------! 332 ELSE ! fixed levels ! (ssh at t-point only) 333 ! !--------------------------! 334 ! 335 IF( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 336 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 337 ! 338 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 339 DO jj = 1, jpj 340 DO ji = 1, jpi ! before <-- now filtered 341 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 342 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 343 END DO 344 END DO 345 ENDIF 346 ! 281 sshn(:,:) = ssha(:,:) ! now <-- after 347 282 ENDIF 348 283 ! … … 354 289 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 355 290 ! 356 END SUBROUTINE ssh_ nxt291 END SUBROUTINE ssh_swp 357 292 358 293 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
