- Timestamp:
- 2010-10-04T15:53:42+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/sshwzv.F90
r2113 r2148 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 7 8 !!---------------------------------------------------------------------- 8 9 … … 27 28 USE diaar5, ONLY : lk_diaar5 28 29 USE iom 29 USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep ! River runoff 30 USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep ! River runoff 30 31 31 32 IMPLICIT NONE … … 38 39 # include "domzgr_substitute.h90" 39 40 # include "vectopt_loop_substitute.h90" 40 41 !!---------------------------------------------------------------------- 42 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 41 !!---------------------------------------------------------------------- 42 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 43 43 !! $Id$ 44 44 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 54 54 !! and update the now vertical coordinate (lk_vvl=T). 55 55 !! 56 !! ** Method : - 57 !! 58 !! - Using the incompressibility hypothesis, the vertical 56 !! ** Method : - Using the incompressibility hypothesis, the vertical 59 57 !! velocity is computed by integrating the horizontal divergence 60 58 !! from the bottom to the surface minus the scale factor evolution. … … 63 61 !! ** action : ssha : after sea surface height 64 62 !! wn : now vertical velocity 65 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 66 !! at u-, v-, f-point s 67 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 63 !! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T) 64 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 65 !! 66 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 68 67 !!---------------------------------------------------------------------- 69 68 USE oce, ONLY : z3d => ta ! use ta as 3D workspace … … 73 72 INTEGER :: ji, jj, jk ! dummy loop indices 74 73 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 75 REAL(wp) :: z2dt, z raur! temporary scalars74 REAL(wp) :: z2dt, z1_2dt, z1_rau0 ! temporary scalars 76 75 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 77 76 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 79 78 80 79 IF( kt == nit000 ) THEN 80 ! 81 81 IF(lwp) WRITE(numout,*) 82 82 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' … … 95 95 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 96 96 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 97 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) &98 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) )99 97 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 100 98 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 101 99 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 102 100 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 103 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) &104 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) )105 101 END DO 106 102 END DO 107 103 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 108 104 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 109 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 105 DO jj = 1, jpjm1 106 DO ji = 1, jpim1 ! NO Vector Opt. 107 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 108 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 109 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 110 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 111 END DO 112 END DO 113 CALL lbc_lnk( sshf_n, 'F', 1. ) 110 114 ENDIF 111 115 ! 112 116 ENDIF 113 117 114 ! !------------------------------ !115 IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)116 ! !------------------------------!118 ! !------------------------------------------! 119 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case) 120 ! !------------------------------------------! 117 121 DO jk = 1, jpkm1 118 fsdept(:,:,jk) = fsdept_n(:,:,jk) 122 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays 119 123 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 120 124 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 121 125 ! 122 fse3t (:,:,jk) = fse3t_n (:,:,jk) 126 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays 123 127 fse3u (:,:,jk) = fse3u_n (:,:,jk) 124 128 fse3v (:,:,jk) = fse3v_n (:,:,jk) … … 128 132 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 129 133 END DO 130 ! ! now ocean depth (at u- and v-points)131 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 134 ! 135 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points) 132 136 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 133 ! 137 ! ! now masked inverse of the ocean depth (at u- and v-points) 134 138 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 135 139 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 136 140 ! 137 ! 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 rnf_dep(ji,jj) = 0. 141 DO jk = 1, rnf_mod_dep(ji,jj) ! recalculates rnf_dep to be the depth 142 rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk) ! in metres to the bottom of the relevant grid box 143 ENDDO 144 ENDDO 145 ENDDO 146 ! 147 ENDIF 148 149 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 150 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 151 152 ! set time step size (Euler/Leapfrog) 153 z2dt = 2. * rdt 141 ENDIF 142 ! 143 144 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 145 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 146 147 z2dt = 2. * rdt ! set time step size (Euler/Leapfrog) 154 148 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 155 156 zraur = 1. / rau0157 149 158 150 ! !------------------------------! … … 163 155 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 164 156 END DO 165 166 157 ! ! Sea surface elevation time stepping 167 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * ( emp(:,:)-rnf(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 158 ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 159 ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b -2*rnf ) = emp 160 z1_rau0 = 0.5 / rau0 161 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) - 2 * rnf(:,:) ) + zhdiv(:,:) ) ) & 162 & * tmask(:,:,1) 168 163 169 164 #if defined key_obc 170 IF 165 IF( Agrif_Root() ) THEN 171 166 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 172 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm)167 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 173 168 ENDIF 174 169 #endif 175 176 170 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 177 171 IF( lk_vvl ) THEN ! (required only in key_vvl case) … … 184 178 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 185 179 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 186 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) &187 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) &188 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) )189 180 END DO 190 181 END DO 191 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 182 ! Boundaries conditions 183 CALL lbc_lnk( sshu_a, 'U', 1. ) 192 184 CALL lbc_lnk( sshv_a, 'V', 1. ) 193 CALL lbc_lnk( sshf_a, 'F', 1. ) 194 ENDIF 195 185 ENDIF 196 186 ! !------------------------------! 197 187 ! ! Now Vertical Velocity ! 198 188 ! !------------------------------! 199 ! ! integrate from the bottom the hor. divergence 200 DO jk = jpkm1, 1, -1 201 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 202 & - ( fse3t_a(:,:,jk) & 203 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 189 z1_2dt = 1.e0 / z2dt 190 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 191 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 192 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 193 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 194 & * tmask(:,:,jk) * z1_2dt 204 195 END DO 205 ! 196 197 ! !------------------------------! 198 ! ! outputs ! 199 ! !------------------------------! 206 200 CALL iom_put( "woce", wn ) ! vertical velocity 207 201 CALL iom_put( "ssh" , sshn ) ! sea surface height 208 202 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 209 IF( lk_diaar5 ) THEN 203 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 204 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 210 205 z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 211 206 DO jk = 1, jpk 212 207 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 213 208 END DO 214 CALL iom_put( "w_masstr" , z3d ) ! vertical mass transport 215 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) ! square of vertical mass transport 216 ENDIF 209 CALL iom_put( "w_masstr" , z3d ) 210 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 211 ENDIF 212 ! 213 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 217 214 ! 218 215 END SUBROUTINE ssh_wzv … … 227 224 !! ssha already computed in ssh_wzv 228 225 !! 229 !! ** Method : - apply Asselin time fiter to now ssh and swap : 230 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 231 !! sshn = ssha 226 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 227 !! from the filter, see Leclair and Madec 2010) and swap : 228 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 229 !! - atfp * rdt * ( emp_b - emp ) / rau0 230 !! sshn = ssha 232 231 !! 233 232 !! ** action : - sshb, sshn : before & now sea surface height 234 233 !! ready for the next time step 235 !!---------------------------------------------------------------------- 236 INTEGER, INTENT( in ) :: kt ! ocean time-step index 237 !! 238 INTEGER :: ji, jj ! dummy loop indices 234 !! 235 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 236 !!---------------------------------------------------------------------- 237 INTEGER, INTENT(in) :: kt ! ocean time-step index 238 !! 239 INTEGER :: ji, jj ! dummy loop indices 240 REAL(wp) :: zec ! temporary scalar 239 241 !!---------------------------------------------------------------------- 240 242 … … 245 247 ENDIF 246 248 247 ! Time filter and swap of the ssh 248 ! ------------------------------- 249 250 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 251 ! ! ---------------------- ! 249 ! !--------------------------! 250 IF( lk_vvl ) THEN ! Variable volume levels ! 251 ! !--------------------------! 252 ! 253 ! ssh at t-, u-, v, f-points 254 !=========================== 252 255 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 253 256 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 254 257 sshu_n(:,:) = sshu_a(:,:) 255 258 sshv_n(:,:) = sshv_a(:,:) 256 sshf_n(:,:) = sshf_a(:,:) 259 DO jj = 1, jpjm1 260 DO ji = 1, jpim1 ! NO Vector Opt. 261 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 262 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 263 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 264 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 265 END DO 266 END DO 267 ! Boundaries conditions 268 CALL lbc_lnk( sshf_n, 'F', 1. ) 257 269 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 270 zec = atfp * rdt / rau0 258 271 DO jj = 1, jpj 259 272 DO ji = 1, jpi ! before <-- now filtered 260 sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 261 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 262 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 263 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 273 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 274 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 264 275 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 265 276 sshu_n(ji,jj) = sshu_a(ji,jj) 266 277 sshv_n(ji,jj) = sshv_a(ji,jj) 267 sshf_n(ji,jj) = sshf_a(ji,jj) 268 END DO 269 END DO 278 END DO 279 END DO 280 DO jj = 1, jpjm1 281 DO ji = 1, jpim1 ! NO Vector Opt. 282 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 283 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 284 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 285 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 286 END DO 287 END DO 288 ! Boundaries conditions 289 CALL lbc_lnk( sshf_n, 'F', 1. ) 290 DO jj = 1, jpjm1 291 DO ji = 1, jpim1 ! NO Vector Opt. 292 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 293 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 294 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 295 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 296 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 297 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 298 END DO 299 END DO 300 ! Boundaries conditions 301 CALL lbc_lnk( sshu_b, 'U', 1. ) 302 CALL lbc_lnk( sshv_b, 'V', 1. ) 270 303 ENDIF 271 ! 272 ELSE ! fixed levels : ssh at t-point only 273 ! ! ------------ ! 304 ! !--------------------------! 305 ELSE ! fixed levels ! 306 ! !--------------------------! 307 ! 308 ! ssh at t-point only 309 !==================== 274 310 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 275 311 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) … … 278 314 DO jj = 1, jpj 279 315 DO ji = 1, jpi ! before <-- now filtered 280 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 316 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 281 317 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 282 318 END DO … … 286 322 ENDIF 287 323 ! 288 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )324 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 289 325 ! 290 326 END SUBROUTINE ssh_nxt
Note: See TracChangeset
for help on using the changeset viewer.
