Changeset 1870 for branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90
- Timestamp:
- 2010-05-12T17:36:00+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1837_MLF/NEMO/OPA_SRC/DYN/sshwzv.F90
r1792 r1870 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 … … 37 38 # include "domzgr_substitute.h90" 38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 40 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 42 42 !! $Id$ 43 43 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 53 53 !! and update the now vertical coordinate (lk_vvl=T). 54 54 !! 55 !! ** Method : - 56 !! 57 !! - Using the incompressibility hypothesis, the vertical 55 !! ** Method : - Using the incompressibility hypothesis, the vertical 58 56 !! velocity is computed by integrating the horizontal divergence 59 57 !! from the bottom to the surface minus the scale factor evolution. … … 62 60 !! ** action : ssha : after sea surface height 63 61 !! wn : now vertical velocity 64 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 65 !! at u-, v-, f-point s 66 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 62 !! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T) 63 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 64 !! 65 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 67 66 !!---------------------------------------------------------------------- 68 67 USE oce, ONLY : z3d => ta ! use ta as 3D workspace … … 78 77 79 78 IF( kt == nit000 ) THEN 79 ! 80 80 IF(lwp) WRITE(numout,*) 81 81 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' … … 111 111 ENDIF 112 112 113 ! !------------------------------ !114 IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)115 ! !------------------------------!113 ! !------------------------------------------! 114 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case) 115 ! !------------------------------------------! 116 116 DO jk = 1, jpkm1 117 fsdept(:,:,jk) = fsdept_n(:,:,jk) 117 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays 118 118 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 119 119 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 120 120 ! 121 fse3t (:,:,jk) = fse3t_n (:,:,jk) 121 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays 122 122 fse3u (:,:,jk) = fse3u_n (:,:,jk) 123 123 fse3v (:,:,jk) = fse3v_n (:,:,jk) … … 127 127 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 128 128 END DO 129 ! ! now ocean depth (at u- and v-points)130 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 129 ! 130 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points) 131 131 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 132 ! 132 ! ! now masked inverse of the ocean depth (at u- and v-points) 133 133 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 134 134 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) … … 136 136 ENDIF 137 137 138 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 139 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 140 141 ! set time step size (Euler/Leapfrog) 142 z2dt = 2. * rdt 138 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 139 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 140 141 z2dt = 2. * rdt ! set time step size (Euler/Leapfrog) 143 142 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 144 143 145 zraur = 1. / rau0146 147 144 ! !------------------------------! 148 145 ! ! After Sea Surface Height ! 149 146 ! !------------------------------! 147 zraur = 0.5 / rau0 148 ! 150 149 zhdiv(:,:) = 0.e0 151 150 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 152 151 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 153 152 END DO 154 155 153 ! ! Sea surface elevation time stepping 156 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)154 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 157 155 158 156 #if defined key_obc 159 IF 157 IF( Agrif_Root() ) THEN 160 158 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 161 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm)159 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 162 160 ENDIF 163 161 #endif 164 165 162 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 166 163 IF( lk_vvl ) THEN ! (required only in key_vvl case) … … 186 183 ! ! Now Vertical Velocity ! 187 184 ! !------------------------------! 188 ! ! integrate from the bottom the hor. divergence 189 DO jk = jpkm1, 1, -1 190 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 191 & - ( fse3t_a(:,:,jk) & 192 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 185 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 186 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 187 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 193 188 END DO 194 ! 189 190 ! !------------------------------! 191 ! ! outputs ! 192 ! !------------------------------! 195 193 CALL iom_put( "woce", wn ) ! vertical velocity 196 194 CALL iom_put( "ssh" , sshn ) ! sea surface height 197 195 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 198 IF( lk_diaar5 ) THEN 196 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 199 197 z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 200 198 DO jk = 1, jpk 201 199 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 202 200 END DO 203 CALL iom_put( "w_masstr" , z3d ) ! vertical mass transport 204 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) ! square of vertical mass transport 205 ENDIF 201 CALL iom_put( "w_masstr" , z3d ) 202 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 203 ENDIF 204 ! 205 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 206 206 ! 207 207 END SUBROUTINE ssh_wzv … … 216 216 !! ssha already computed in ssh_wzv 217 217 !! 218 !! ** Method : - apply Asselin time fiter to now ssh and swap : 219 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 220 !! sshn = ssha 218 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 219 !! from the filter, see Leclair and Madec 2010) and swap : 220 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 221 !! - atfp * rdt * ( emp_b - emp ) / rau0 222 !! sshn = ssha 221 223 !! 222 224 !! ** action : - sshb, sshn : before & now sea surface height 223 225 !! ready for the next time step 224 !!---------------------------------------------------------------------- 225 INTEGER, INTENT( in ) :: kt ! ocean time-step index 226 !! 227 INTEGER :: ji, jj ! dummy loop indices 226 !! 227 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 228 !!---------------------------------------------------------------------- 229 INTEGER, INTENT(in) :: kt ! ocean time-step index 230 !! 231 INTEGER :: ji, jj ! dummy loop indices 232 REAL(wp) :: zec ! temporary scalare 228 233 !!---------------------------------------------------------------------- 229 234 … … 234 239 ENDIF 235 240 236 ! Time filter and swap of the ssh 237 ! ------------------------------- 238 239 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 240 ! ! ---------------------- ! 241 ! !--------------------------! 242 IF( lk_vvl ) THEN ! Variable volume levels ! ssh at t-, u-, v, f-points 243 ! !--------------------------! 241 244 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 242 245 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) … … 247 250 DO jj = 1, jpj 248 251 DO ji = 1, jpi ! before <-- now filtered 249 sshb (ji,jj) = sshn (ji,jj)+ atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) )252 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 250 253 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 251 254 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) … … 258 261 END DO 259 262 ENDIF 260 ! 261 ELSE ! fixed levels :ssh at t-point only262 ! ! ------------!263 ! !--------------------------! 264 ELSE ! fixed levels ! ssh at t-point only 265 ! !--------------------------! 263 266 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 264 267 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 265 268 ! 266 269 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 270 zec = atfp * rdt / rau0 267 271 DO jj = 1, jpj 268 272 DO ji = 1, jpi ! before <-- now filtered 269 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(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) 270 275 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 271 276 END DO … … 274 279 ! 275 280 ENDIF 276 ! 277 IF(ln_ctl) CALL prt_ctl(tab2d_1=sshb , clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 281 282 ! time filter with global conservation correction and array swap 283 sshbb(:,:) = sshb(:,:) 284 sshb (:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + zssha(:,:) ) & 285 & - zfact * 286 sshn (:,:) = zssha(:,:) 287 empb (:,:) = emp(:,:) 288 289 290 ! 291 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 278 292 ! 279 293 END SUBROUTINE ssh_nxt
Note: See TracChangeset
for help on using the changeset viewer.
