Changeset 1438 for trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
- Timestamp:
- 2009-05-11T16:34:47+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/wzvmod.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
r1241 r1438 1 1 MODULE wzvmod 2 !! MODULE sshwzv 2 3 !!============================================================================== 3 !! *** MODULE wzvmod***4 !! Ocean d iagnostic variable :vertical velocity4 !! *** MODULE sshwzv *** 5 !! Ocean dynamics : sea surface height and vertical velocity 5 6 !!============================================================================== 6 !! History : 5.0 ! 90-10 (C. Levy, G. Madec) Original code7 !! 7.0 ! 96-01 (G. Madec) Statement function for e38 !! 8.5 ! 02-07 (G. Madec) Free form, F90 9 !! " ! 07-07 (D. Storkey) Zero zhdiv at open boundary (BDY)10 !! ----------------------------------------------------------------------11 !! wzv : Compute the vertical velocity12 !! ----------------------------------------------------------------------13 !! * Modules used7 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 8 !!---------------------------------------------------------------------- 9 10 !!---------------------------------------------------------------------- 11 !! ssh_wzv : after ssh & now vertical velocity 12 !! ssh_nxt : filter ans swap the ssh arrays 13 !! ssh_rst : read/write ssh restart fields in the ocean restart file 14 !!---------------------------------------------------------------------- 14 15 USE oce ! ocean dynamics and tracers variables 15 16 USE dom_oce ! ocean space and time domain variables 16 17 USE sbc_oce ! surface boundary condition: ocean 17 18 USE domvvl ! Variable volume 19 USE iom ! I/O library 20 USE restart ! only for lrst_oce 18 21 USE in_out_manager ! I/O manager 19 22 USE prtctl ! Print control 20 23 USE phycst 21 USE bdy_oce ! unstructured open boundaries22 24 USE lbclnk ! ocean lateral boundary condition (or mpp link) 23 25 USE obc_par ! open boundary cond. parameter … … 27 29 PRIVATE 28 30 29 !! * Routine accessibility30 PUBLIC wzv ! routine called by step.F90 and inidtr.F9031 PUBLIC ssh_wzv ! called by step.F90 32 PUBLIC ssh_nxt ! called by step.F90 31 33 32 34 !! * Substitutions 33 35 # include "domzgr_substitute.h90" 34 !!---------------------------------------------------------------------- 35 !! OPA 9.0 , LOCEAN-IPSL (2005) 36 # include "vectopt_loop_substitute.h90" 37 38 !!---------------------------------------------------------------------- 39 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 36 40 !! $Id$ 37 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 40 44 CONTAINS 41 45 42 SUBROUTINE wzv( kt )43 !!---------------------------------------------------------------------- 44 !! *** ROUTINEwzv ***45 !! 46 !! ** Purpose : Compute the now vertical velocity after the array swap47 !! 48 !! ** Method : Using the incompressibility hypothesis, the vertical49 !! velocity is computed by integrating the horizontal divergence50 !! from the bottom to the surface.51 !! The boundary conditions are w=0 at the bottom (no flux) and,52 !! in regid-lid case, w=0 at the sea surface.53 !! 54 !! ** action : wn array : the now vertical velocity55 !! ----------------------------------------------------------------------56 !! * Arguments57 INTEGER, INTENT( in ) :: kt ! ocean time-step index58 59 !! * Local declarations60 INTEGER :: jk ! dummy loop indices61 !! Variable volume62 INTEGER :: ji, jj ! dummy loop indices63 REAL(wp) :: z2dt, zraur ! temporary scalar64 REAL(wp), DIMENSION (jpi,jpj) :: zssha, zun, zvn, zhdiv65 #if defined key_bdy 66 INTEGER :: jgrd, jb! temporary scalars67 #endif 46 SUBROUTINE ssh_wzv( kt ) 47 !!---------------------------------------------------------------------- 48 !! *** ROUTINE ssh_wzv *** 49 !! 50 !! ** Purpose : compute the after ssh (ssha), the now vertical velocity 51 !! and update the now vertical coordinate (lk_vvl=T). 52 !! 53 !! ** Method : - 54 !! 55 !! - Using the incompressibility hypothesis, the vertical 56 !! velocity is computed by integrating the horizontal divergence 57 !! from the bottom to the surface minus the scale factor evolution. 58 !! The boundary conditions are w=0 at the bottom (no flux) and. 59 !! 60 !! ** action : ssha : after sea surface height 61 !! wn : now vertical velocity 62 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 63 !! at u-, v-, f-point s 64 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 65 !!---------------------------------------------------------------------- 66 INTEGER, INTENT(in) :: kt ! time step 67 !! 68 INTEGER :: ji, jj, jk ! dummy loop indices 69 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 70 REAL(wp) :: z2dt, zraur ! temporary scalars 71 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 68 72 !!---------------------------------------------------------------------- 69 73 70 74 IF( kt == nit000 ) THEN 71 75 IF(lwp) WRITE(numout,*) 72 IF(lwp) WRITE(numout,*) 'wzv : vertical velocity from continuity eq.' 73 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 74 75 ! bottom boundary condition: w=0 (set once for all) 76 wn(:,:,jpk) = 0.e0 77 ENDIF 78 79 IF( lk_vvl ) THEN ! Variable volume 80 ! 81 z2dt = 2. * rdt ! time step: leap-frog 82 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 83 zraur = 1. / rauw 84 85 ! Vertically integrated quantities 86 ! -------------------------------- 87 zun(:,:) = 0.e0 88 zvn(:,:) = 0.e0 89 ! 90 DO jk = 1, jpkm1 ! Vertically integrated transports (now) 91 zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 92 zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 93 END DO 94 95 ! Horizontal divergence of barotropic transports 96 !-------------------------------------------------- 97 zhdiv(:,:) = 0.e0 98 DO jj = 2, jpjm1 99 DO ji = 2, jpim1 ! vector opt. 100 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun(ji ,jj ) & 101 & - e2u(ji-1,jj ) * zun(ji-1,jj ) & 102 & + e1v(ji ,jj ) * zvn(ji ,jj ) & 103 & - e1v(ji ,jj-1) * zvn(ji ,jj-1) ) & 104 & / ( e1t(ji,jj) * e2t(ji,jj) ) 76 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' 77 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 78 ! 79 CALL ssh_rst( nit000, 'READ' ) ! read or initialize sshb and sshn 80 ! 81 wn(:,:,jpk) = 0.e0 ! bottom boundary condition: w=0 (set once for all) 82 ! 83 IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only) 84 DO jj = 1, jpjm1 85 DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible 86 zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) 87 zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) 88 zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) 89 sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 90 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 91 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 92 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 93 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) & 94 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) 95 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 96 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 97 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 98 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 99 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) & 100 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) 101 END DO 102 END DO 103 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 104 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 105 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 106 ENDIF 107 ! 108 ENDIF 109 110 ! set time step size (Euler/Leapfrog) 111 z2dt = 2. * rdt 112 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 113 114 zraur = 1. / rauw 115 116 ! !------------------------------! 117 ! ! After Sea Surface Height ! 118 ! !------------------------------! 119 zhdiv(:,:) = 0.e0 120 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 121 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 122 END DO 123 124 ! ! Sea surface elevation time stepping 125 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 126 127 #if defined key_obc 128 # if defined key_agrif 129 IF ( Agrif_Root() ) THEN 130 # endif 131 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 132 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 133 # if defined key_agrif 134 ENDIF 135 # endif 136 #endif 137 138 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 139 IF( lk_vvl ) THEN ! (required only in key_vvl case) 140 DO jj = 1, jpjm1 141 DO ji = 1, fs_jpim1 ! Vector Opt. 142 sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 143 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 144 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 145 sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 146 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 147 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 148 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) & ! Caution : fmask not used 149 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) & 150 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) ) 105 151 END DO 106 152 END DO 107 108 #if defined key_obc && ( defined key_dynspg_exp || defined key_dynspg_ts ) 109 ! open boundaries (div must be zero behind the open boundary) 110 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 111 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east 112 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west 113 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 114 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south 115 #endif 116 117 #if defined key_bdy 118 jgrd=1 !: tracer grid. 119 DO jb = 1, nblenrim(jgrd) 120 ji = nbi(jb,jgrd) 121 jj = nbj(jb,jgrd) 122 zhdiv(ji,jj) = 0.e0 153 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 154 CALL lbc_lnk( sshv_a, 'V', 1. ) 155 CALL lbc_lnk( sshf_a, 'F', 1. ) 156 ENDIF 157 158 ! !------------------------------! 159 ! ! Now Vertical Velocity ! 160 ! !------------------------------! 161 ! ! integrate from the bottom the hor. divergence 162 DO jk = jpkm1, 1, -1 163 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 164 & - ( fse3t_a(:,:,jk) & 165 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 166 END DO 167 168 ! !------------------------------! 169 ! ! Update Now Vertical coord. ! 170 ! !------------------------------! 171 IF( lk_vvl ) THEN ! only in vvl case) 172 ! ! now local depth and scale factors (stored in fse3. arrays) 173 DO jk = 1, jpkm1 174 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! depths 175 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 176 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 177 ! 178 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors 179 fse3u (:,:,jk) = fse3u_n (:,:,jk) 180 fse3v (:,:,jk) = fse3v_n (:,:,jk) 181 fse3f (:,:,jk) = fse3f_n (:,:,jk) 182 fse3w (:,:,jk) = fse3w_n (:,:,jk) 183 fse3uw(:,:,jk) = fse3uw_n(:,:,jk) 184 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 123 185 END DO 124 #endif 125 126 CALL lbc_lnk( zhdiv, 'T', 1. ) 127 128 ! Sea surface elevation time stepping 129 ! ----------------------------------- 130 zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 131 132 ! Vertical velocity computed from bottom 133 ! -------------------------------------- 134 DO jk = jpkm1, 1, -1 135 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) & 136 & - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt 137 END DO 138 139 ELSE ! Fixed volume 140 141 ! Vertical velocity computed from bottom 142 ! -------------------------------------- 143 DO jk = jpkm1, 1, -1 144 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) 145 END DO 146 147 ENDIF 148 149 IF(ln_ctl) CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 - : ', mask1=wn) 150 151 END SUBROUTINE wzv 186 ! ! ocean depth (at u- and v-points) 187 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 188 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 189 ! ! masked inverse of the ocean depth (at u- and v-points) 190 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 191 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 192 193 ENDIF 194 ! 195 END SUBROUTINE ssh_wzv 196 197 198 SUBROUTINE ssh_nxt( kt ) 199 !!---------------------------------------------------------------------- 200 !! *** ROUTINE ssh_nxt *** 201 !! 202 !! ** Purpose : achieve the sea surface height time stepping by 203 !! applying Asselin time filter and swapping the arrays 204 !! ssha already computed in ssh_wzv 205 !! 206 !! ** Method : - apply Asselin time fiter to now ssh and swap : 207 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 208 !! sshn = ssha 209 !! 210 !! ** action : - sshb, sshn : before & now sea surface height 211 !! ready for the next time step 212 !!---------------------------------------------------------------------- 213 INTEGER, INTENT( in ) :: kt ! ocean time-step index 214 !! 215 INTEGER :: ji, jj ! dummy loop indices 216 !!---------------------------------------------------------------------- 217 218 IF( kt == nit000 ) THEN 219 IF(lwp) WRITE(numout,*) 220 IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' 221 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 222 ENDIF 223 224 ! Time filter and swap of the ssh 225 ! ------------------------------- 226 227 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 228 ! ! ---------------------- ! 229 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 230 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 231 sshu_n(:,:) = sshu_a(:,:) 232 sshv_n(:,:) = sshv_a(:,:) 233 sshf_n(:,:) = sshf_a(:,:) 234 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 235 DO jj = 1, jpj 236 DO ji = 1, jpi ! before <-- now filtered 237 sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 238 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 239 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 240 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 241 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 242 sshu_n(ji,jj) = sshu_a(ji,jj) 243 sshv_n(ji,jj) = sshv_a(ji,jj) 244 sshf_n(ji,jj) = sshf_a(ji,jj) 245 END DO 246 END DO 247 ENDIF 248 ! 249 ELSE ! fixed levels : ssh at t-point only 250 ! ! ------------ ! 251 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 252 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 253 ! 254 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 255 DO jj = 1, jpj 256 DO ji = 1, jpi ! before <-- now filtered 257 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 258 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 259 END DO 260 END DO 261 ENDIF 262 ! 263 ENDIF 264 265 ! ! write filtered free surface arrays in restart file 266 IF( lrst_oce ) CALL ssh_rst( kt, 'WRITE' ) 267 ! 268 IF(ln_ctl) CALL prt_ctl(tab2d_1=sshb , clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 269 ! 270 END SUBROUTINE ssh_nxt 271 272 273 SUBROUTINE ssh_rst( kt, cdrw ) 274 !!--------------------------------------------------------------------- 275 !! *** ROUTINE ssh_rst *** 276 !! 277 !! ** Purpose : Read or write Sea Surface Height arrays in restart file 278 !! 279 !! ** action : - cdrw = READ : sshb, sshn read in ocean restart file 280 !! - cdrw = WRITE : sshb, sshn written in ocean restart file 281 !!---------------------------------------------------------------------- 282 INTEGER , INTENT(in) :: kt ! ocean time-step 283 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 284 !!---------------------------------------------------------------------- 285 ! 286 IF( TRIM(cdrw) == 'READ' ) THEN 287 IF( iom_varid( numror, 'sshn', ldstop = .FALSE. ) > 0 ) THEN 288 CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) ) 289 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) ) 290 IF( neuler == 0 ) sshb(:,:) = sshn(:,:) 291 ELSE 292 IF( nn_rstssh == 1 ) THEN 293 sshb(:,:) = 0.e0 294 sshn(:,:) = 0.e0 295 ENDIF 296 ENDIF 297 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 298 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb(:,:) ) 299 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn(:,:) ) 300 ENDIF 301 ! 302 END SUBROUTINE ssh_rst 152 303 153 304 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
