Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
- Property svn:eol-style deleted
r2486 r2528 5 5 !!============================================================================== 6 6 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 7 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 8 !! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 7 10 !!---------------------------------------------------------------------- 8 11 … … 16 19 USE domvvl ! Variable volume 17 20 USE divcur ! hor. divergence and curl (div & cur routines) 18 USE cla_div ! cross land: hor. divergence (div_cla routine)19 21 USE iom ! I/O library 20 22 USE restart ! only for lrst_oce … … 25 27 USE obc_par ! open boundary cond. parameter 26 28 USE obc_oce 29 USE bdy_oce 27 30 USE diaar5, ONLY : lk_diaar5 28 31 USE iom 32 USE sbcrnf, ONLY : h_rnf, nk_rnf ! River runoff 33 #if defined key_agrif 34 USE agrif_opa_update 29 35 USE agrif_opa_interp 30 USE agrif_opa_update 36 #endif 37 #if defined key_asminc 38 USE asminc ! Assimilation increment 39 #endif 31 40 32 41 IMPLICIT NONE … … 39 48 # include "domzgr_substitute.h90" 40 49 # include "vectopt_loop_substitute.h90" 41 42 !!---------------------------------------------------------------------- 43 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 50 !!---------------------------------------------------------------------- 51 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 44 52 !! $Id$ 45 !! Software governed by the CeCILL licence ( modipsl/doc/NEMO_CeCILL.txt)53 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 46 54 !!---------------------------------------------------------------------- 47 55 … … 55 63 !! and update the now vertical coordinate (lk_vvl=T). 56 64 !! 57 !! ** Method : - 58 !! 59 !! - Using the incompressibility hypothesis, the vertical 65 !! ** Method : - Using the incompressibility hypothesis, the vertical 60 66 !! velocity is computed by integrating the horizontal divergence 61 67 !! from the bottom to the surface minus the scale factor evolution. … … 64 70 !! ** action : ssha : after sea surface height 65 71 !! wn : now vertical velocity 66 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 67 !! at u-, v-, f-point s 68 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 72 !! 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-points 74 !! 75 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 69 76 !!---------------------------------------------------------------------- 70 77 USE oce, ONLY : z3d => ta ! use ta as 3D workspace … … 74 81 INTEGER :: ji, jj, jk ! dummy loop indices 75 82 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 76 REAL(wp) :: z2dt, z raur! temporary scalars83 REAL(wp) :: z2dt, z1_2dt, z1_rau0 ! temporary scalars 77 84 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 78 85 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 80 87 81 88 IF( kt == nit000 ) THEN 89 ! 82 90 IF(lwp) WRITE(numout,*) 83 91 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' … … 96 104 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 97 105 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 98 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) &99 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) )100 106 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 101 107 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 102 108 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 103 109 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 104 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) &105 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) )106 110 END DO 107 111 END DO 108 112 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 109 113 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 110 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 114 DO jj = 1, jpjm1 115 DO ji = 1, jpim1 ! NO Vector Opt. 116 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 117 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 118 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 119 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 120 END DO 121 END DO 122 CALL lbc_lnk( sshf_n, 'F', 1. ) 111 123 ENDIF 112 124 ! 113 125 ENDIF 114 126 115 ! !------------------------------ !116 IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)117 ! !------------------------------!127 ! !------------------------------------------! 128 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case) 129 ! !------------------------------------------! 118 130 DO jk = 1, jpkm1 119 fsdept(:,:,jk) = fsdept_n(:,:,jk) 131 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays 120 132 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 121 133 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 122 134 ! 123 fse3t (:,:,jk) = fse3t_n (:,:,jk) 135 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays 124 136 fse3u (:,:,jk) = fse3u_n (:,:,jk) 125 137 fse3v (:,:,jk) = fse3v_n (:,:,jk) … … 129 141 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 130 142 END DO 131 ! ! now ocean depth (at u- and v-points)132 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 143 ! 144 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points) 133 145 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 134 ! 146 ! ! now masked inverse of the ocean depth (at u- and v-points) 135 147 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 136 148 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 137 ! 138 ENDIF 139 140 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 141 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 142 143 ! set time step size (Euler/Leapfrog) 144 z2dt = 2. * rdt 149 ! 150 ENDIF 151 ! 152 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 153 ! 154 z2dt = 2. * rdt ! set time step size (Euler/Leapfrog) 145 155 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 146 147 zraur = 1. / rau0148 156 149 157 ! !------------------------------! … … 154 162 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 155 163 END DO 156 157 164 ! ! Sea surface elevation time stepping 158 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 165 ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 166 ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp 167 z1_rau0 = 0.5 / rau0 168 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) & 169 & * tmask(:,:,1) 159 170 160 171 #if defined key_agrif … … 162 173 #endif 163 174 #if defined key_obc 164 IF 175 IF( Agrif_Root() ) THEN 165 176 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 166 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 167 ENDIF 177 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 178 ENDIF 179 #endif 180 #if defined key_bdy 181 ssha(:,:) = ssha(:,:) * bdytmask(:,:) 182 CALL lbc_lnk( ssha, 'T', 1. ) 168 183 #endif 169 184 … … 178 193 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 179 194 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 180 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) &181 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) &182 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) )183 195 END DO 184 196 END DO 185 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 197 ! Boundaries conditions 198 CALL lbc_lnk( sshu_a, 'U', 1. ) 186 199 CALL lbc_lnk( sshv_a, 'V', 1. ) 187 CALL lbc_lnk( sshf_a, 'F', 1. ) 188 ENDIF 200 ENDIF 201 ! Include the IAU weighted SSH increment 202 #if defined key_asminc 203 IF( ( lk_asminc ).AND.( ln_sshinc ).AND.( ln_asmiau ) ) THEN 204 CALL ssh_asm_inc( kt ) 205 ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 206 ENDIF 207 #endif 189 208 190 209 ! !------------------------------! 191 210 ! ! Now Vertical Velocity ! 192 211 ! !------------------------------! 193 ! ! integrate from the bottom the hor. divergence 194 DO jk = jpkm1, 1, -1 195 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 196 & - ( fse3t_a(:,:,jk) & 197 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 212 z1_2dt = 1.e0 / z2dt 213 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 214 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 215 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 216 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 217 & * tmask(:,:,jk) * z1_2dt 218 #if defined key_bdy 219 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) 220 #endif 198 221 END DO 199 ! 222 223 ! !------------------------------! 224 ! ! outputs ! 225 ! !------------------------------! 200 226 CALL iom_put( "woce", wn ) ! vertical velocity 201 227 CALL iom_put( "ssh" , sshn ) ! sea surface height 202 228 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 203 IF( lk_diaar5 ) THEN 229 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 230 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 204 231 z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 205 232 DO jk = 1, jpk 206 233 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 207 234 END DO 208 CALL iom_put( "w_masstr" , z3d ) ! vertical mass transport 209 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) ! square of vertical mass transport 210 ENDIF 235 CALL iom_put( "w_masstr" , z3d ) 236 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 237 ENDIF 238 ! 239 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 211 240 ! 212 241 END SUBROUTINE ssh_wzv … … 221 250 !! ssha already computed in ssh_wzv 222 251 !! 223 !! ** Method : - apply Asselin time fiter to now ssh and swap : 224 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 225 !! sshn = ssha 252 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 253 !! from the filter, see Leclair and Madec 2010) and swap : 254 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 255 !! - atfp * rdt * ( emp_b - emp ) / rau0 256 !! sshn = ssha 226 257 !! 227 258 !! ** action : - sshb, sshn : before & now sea surface height 228 259 !! ready for the next time step 229 !!---------------------------------------------------------------------- 230 INTEGER, INTENT( in ) :: kt ! ocean time-step index 231 !! 232 INTEGER :: ji, jj ! dummy loop indices 260 !! 261 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 262 !!---------------------------------------------------------------------- 263 INTEGER, INTENT(in) :: kt ! ocean time-step index 264 !! 265 INTEGER :: ji, jj ! dummy loop indices 266 REAL(wp) :: zec ! temporary scalar 233 267 !!---------------------------------------------------------------------- 234 268 … … 239 273 ENDIF 240 274 241 ! Time filter and swap of the ssh 242 ! ------------------------------- 243 244 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 245 ! ! ---------------------- ! 275 ! !--------------------------! 276 IF( lk_vvl ) THEN ! Variable volume levels ! 277 ! !--------------------------! 278 ! 279 ! ssh at t-, u-, v, f-points 280 !=========================== 246 281 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 247 282 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 248 283 sshu_n(:,:) = sshu_a(:,:) 249 284 sshv_n(:,:) = sshv_a(:,:) 250 sshf_n(:,:) = sshf_a(:,:) 285 DO jj = 1, jpjm1 286 DO ji = 1, jpim1 ! NO Vector Opt. 287 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 288 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 289 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 290 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 291 END DO 292 END DO 293 ! Boundaries conditions 294 CALL lbc_lnk( sshf_n, 'F', 1. ) 251 295 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 296 zec = atfp * rdt / rau0 252 297 DO jj = 1, jpj 253 298 DO ji = 1, jpi ! before <-- now filtered 254 sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 255 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 256 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 257 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 299 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 300 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 258 301 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 259 302 sshu_n(ji,jj) = sshu_a(ji,jj) 260 303 sshv_n(ji,jj) = sshv_a(ji,jj) 261 sshf_n(ji,jj) = sshf_a(ji,jj) 262 END DO 263 END DO 304 END DO 305 END DO 306 DO jj = 1, jpjm1 307 DO ji = 1, jpim1 ! NO Vector Opt. 308 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 309 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 310 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 311 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 312 END DO 313 END DO 314 ! Boundaries conditions 315 CALL lbc_lnk( sshf_n, 'F', 1. ) 316 DO jj = 1, jpjm1 317 DO ji = 1, jpim1 ! NO Vector Opt. 318 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 319 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 320 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 321 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 322 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 323 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 324 END DO 325 END DO 326 ! Boundaries conditions 327 CALL lbc_lnk( sshu_b, 'U', 1. ) 328 CALL lbc_lnk( sshv_b, 'V', 1. ) 264 329 ENDIF 265 ! 266 ELSE ! fixed levels : ssh at t-point only 267 ! ! ------------ ! 330 ! !--------------------------! 331 ELSE ! fixed levels ! 332 ! !--------------------------! 333 ! 334 ! ssh at t-point only 335 !==================== 268 336 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 269 337 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) … … 272 340 DO jj = 1, jpj 273 341 DO ji = 1, jpi ! before <-- now filtered 274 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 342 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 275 343 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 276 344 END DO … … 280 348 ENDIF 281 349 ! 350 ! Update velocity at AGRIF zoom boundaries 282 351 #if defined key_agrif 283 ! Update velocity at AGRIF zoom boundaries 284 IF (.NOT.Agrif_Root()) CALL Agrif_Update_Dyn( kt ) 285 #endif 286 287 IF(ln_ctl) CALL prt_ctl(tab2d_1=sshb , clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 352 IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt ) 353 #endif 354 ! 355 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 288 356 ! 289 357 END SUBROUTINE ssh_nxt
Note: See TracChangeset
for help on using the changeset viewer.
